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ABSTRACT 

In  this  thesis  single-degree-of-freedom  torsional  airfoil  flutter  is  investigated 
using  an  incompressible  potential  flow  code,  a  compressible  inviscid  Euler  code  and  a 
compressible  viscous  Navier-Stokes  code.  It  is  found  that  the  classical  linearized 
incompressible  and  compressible  flow  theories  yield  unconservative  flutter  estimates. 
The  computations  based  on  the  non-linear  codes  show  for  NACA  0006,  NACA  0009, 
NACA  0012  and  NACA  0015  airfoils  that  the  regions  of  torsional  flutter  instability 
increase  as  the  airfoil  thickness  and  the  flight  Mach  number  is  increased.  On  the  other 
hand,  the  comparison  of  the  flutter  boundaries  computed  with  the  viscous  Navier-Stokes 
code  versus  the  inviscid  Euler  code  shows  that  the  effect  of  viscosity  is  stabilizing.  Also, 
the  computed  flutter  boundaries  display  the  effect  of  pitch  axis  location  on  flutter.  Axis 
locations  in  the  range  between  half  a  chord  upstream  of  the  leading  edge  of  the  airfoil  and 
the  leading  edge  are  most  prone  to  induce  flutter.  Axis  locations  downstream  of  the 
quarter  chord  are  flutter  free. 
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I.  INTRODUCTION 


A.  BACKGROUND 

The  phenomenon  of  airfoil  flutter  was  first  encountered  on  World  War  I  airplanes. 
The  danger  to  flight  safety  posed  by  this  phenomenon  stimulated  the  development  of 
quite  sophisticated  flutter  analysis  and  testing  methods  in  the  ensuing  years.  References 
1,3, 4, 6  and  7  give  good  reviews  of  these  developments,  where  it  is  also  explained  that  the 
flutter  phenomenon  is  caused  by  the  interaction  of  the  aerodynamic,  inertia  and  elastic 
forces.  Until  very  recently,  important  simplifying  assumptions  had  to  be  introduced  in 
order  to  make  the  problem  mathematically  tractable.  Linearized  aerodynamic  analysis 
methods  were  among  the  most  important  simplifications  used  in  the  past.  For  the  analysis 
of  low  speed  flutter  phenomena,  Theodorsen  developed  an  incompressible  flow  theory 
for  oscillating  flat  plates,  which  was  then  extended  to  compressible  subsonic  flow.  This 
incompressible  oscillatory  flat  plate  theory  remains  the  standard  aerodynamic  analysis 
tool  for  flutter  calculations.  A  second  important  simplification  was  also  introduced  by 
Theodorsen  when  he  proposed  to  perform  the  flutter  analysis  in  the  frequency  domain.  In 
this  approach,  the  aerodynamic  forces  need  to  be  computed  only  for  the  special  case 
where  the  airfoil  is  assumed  to  execute  a  purely  harmonic  oscillation  at  constant  small 
amplitude. 

Flutter  analyses  based  on  linearized  aerodynamic  theory  in  the  frequency  domain 
have  serious  shortcomings.  Thin  airfoil  (flat  plate)  theory  makes  it  impossible  to  account 
for  the  effect  of  airfoil  geometry  on  flutter.  Frequency-domain  methods  provide  little 
information  about  the  physics  of  the  flutter  problem  because  the  decay  of  the  airfoil 
motion  below  the  critical  flutter  speed  and  the  divergence  of  the  motion  above  this  speed 
cannot  be  obtained  as  part  of  frequency  domain  flutter  analysis. 

The  recent  rapid  advances  in  computational  fluid  dynamics  (CFD),  however, 
made  it  possible  to  replace  linearized  aerodynamic  methods  by  solutions  based  on  the  full 
non-linear  Euler  or  Navier-Stokes  equations  for  inviscid  or  viscous  compressible  flow. 
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As  an  additional  advantage,  these  modem  CFD  methods  integrate  the  governing  flow 
equations  in  time.  Therefore,  they  can  easily  and  naturally  be  combined  with  the 
equations  of  motion  of  the  airfoil.  This  leads  to  the  following  specific  objective  for  the 
work  attempted  in  this  investigation. 


B.  OBJECTIVE 

A  two-dimensional  compressible  Euler  and  Navier-Stokes  flow  solver  is  coupled 
with  a  one  degree-of-freedom  structural  model  for  the  time  domain  analysis  of  an  airfoil 
which  is  free  to  oscillate  about  a  specific  pivot  point.  The  effect  of  airfoil  thickness  on 
torsional  flutter  is  to  be  determined  as  a  function  of  pivot  location  and  subsonic  Mach 
number.  Also,  viscous  flow  effects  on  flutter  are  to  be  determined  by  comparing  Euler 
and  Navier-Stokes  computed  flutter  boundaries.  Furthermore,  at  low  Mach  numbers  the 
Euler  computed  flutter  boundaries  are  to  be  compared  with  the  flutter  boundaries 
computed  with  the  incompressible  unsteady  panel  code  UPOT. 
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II.  GOVERNING  EQUATIONS 


The  equations  that  describe  the  compressible  viscous  fluid  flow  around  a  body  are 
the  continuity,  the  momentum  and  the  energy  equations.  The  flow  around  the  body  can  be 
computed  by  the  simultaneous  solution  of  these  equations.  The  conservation-law  and  the 
vector  form  of  the  compressible  Reynolds-averaged  Navier-Stokes  equation  is  presented 
in  this  chapter.  A  detailed  derivation  can  be  found  in  [Ref  8]. 


A.  CONTINUITY  EQUATION 


The  continuity  equation  expresses  the  conservation-of-mass  law  applied  to  a 
fluid  passing  through  a  control  volume  fixed  in  space 

^  +  V(pV)  =  0  (2.1) 

dt 

where  p  is  the  fluid  density  and  V  is  the  fluid  velocity.  Eq.(2.1)  states  that  the  net  mass 
flux  through  a  control  volume  bounding  surface  must  be  equal  to  the  time  rate  of  change 
of  the  mass  inside  the  control  volume.  For  two-dimensional  Cartesian  flow  this  equation 
reads 

-^•  +  —  (y0M)+— (/7W)=.°  (2.2) 

where  u  and  w  are  velocity  components  along  the  x  and  z  directions,  respectively. 


B.  MOMENTUM  EQUATION 

The  momentum  equation  expresses  Newton’s  second  law  as  applied  to  a  fluid 
element  flowing  relative  to  a  space-fixed  coordinate  system.  The  x  and  z  direction 
momentum  equations  are: 
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(2.3) 


^l  +  V(pW)  =  ff  +  VYliJ 

at 

The  first  term  in  Eq.(2.3)  represents  the  time  rate  of  change  of  momentum  per  unit 
volume  in  the  control  volume.  The  second  term  represents  the  momentum  net  flux 
through  the  bounding  surface  of  the  control  volume.  Also,  f  is  the  body  force  per  unit 

volume  and  y  is  the  stress  tensor  given  by: 

rL=-f>'W 


dui  2  ^  duk 
dxj  3  ‘J  dx 


k  j 


(2.4) 


where  i,j,k  =  1,2,3  and  is  the  Kronecker  delta. 

By  substituting  Eq.(2.4)  into  (2.3)  for  flow  in  a  two-dimensional  Cartesian 
coordinate  system  we  obtain: 
dpu  dpu  dpu  .  dp  d 


.  ,ndu  3w 

d 

f  dw  du ") 

2/3//(2—  — 
dx  dz  _ 

+  dz 

M 

dx  dz )_ 

(2.5) 


dpw  dpw  dpw  dp  d 

a+u&+war=fif’-*+& 


2/3// 


(  dw  di^ 
dz  dx. 


dx 


dw  ■  du 


These  equations  are  known  as  the  Navier-Stokes  equations. 


C.  ENERGY  EQUATION 

The  energy  equation  is  derived  by  applying  the  first  law  of  thermodynamics  (rate 
of  change  of  energy  =  net  heat  flux  into  particle  +  rate  of  work  done  on  particle). 

p(f  u  +  f  w)  +  ^-(eu  +  pu-urxx-wixz)+  (2.6) 

dt  dt  dx 

—  (ew  +  pw-  wt a  -urxz)-0 
dz 

where  e  is  the  total  energy  per  unit  volume. 
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D.  VECTOR  FORM  OF  FLOW  EQUATIONS 

The  above  equations  can  be  rewritten  in  non-dimensionalized  vector  form  as 


where 


Q  = 


,  9F  ,  3G_  1  3F.  |  8G 
Bt  Bx  Bz  Re  dx  dz 


~  p~ 

pu 

pw 

pu 

F  = 

pu 2  +  p 

G- 

puw 

2 

pw 

puw 

pw  +  p 

e 

(e  +  p)u 

(e  +  p)w 

(2.7) 


and 


F  = 


'O' 

"O' 

T 

XX 

Gv  = 

xz 

T 

T 

xz 

7Z 

A. 

J*. 

Txx  =~M(ux-l/2wz) 


Tzz=^M(^z~l/2ux) 


=M(wx  +uz) 

/4  =  UT^  +  WTC  + 


*  a2 


g4=UTxz+WTzz  + 


Re  = 


Pr(r-l) 

a 

Pr(7-1)‘ 
UL 


is  the  free  stream  velocity  and  L  is  the  reference  length.  The  pressure  is  related 
to  the  other  variables  by 
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In  the  above  equations  the  ratio  of  the  specific  heats,  yis  1.4  and  a  is  the  local  speed  of 
sound.  The  density  is  non-dimensionalized  by  the  free  stream  density  p„ ,  the  velocities 
by  the  free  stream  speed  of  sound  ,  and  the  total  energy  by  pmaj  . 


E.  TURBULENCE  MODEL 

The  unsteady  Navier-Stokes  equations  can  completely  model  the  fluid  flow,  but 
the  computational  calculations  of  turbulent  flows  for  realistic  geometries  at  high 
Reynolds  numbers  demand  very  high  grid  densities  and  very  small  time  steps.  Therefore, 
in  order  to  compute  turbulent  flows  for  configurations  of  practical  interest,  turbulence 
modeling  is  used.  Turbulence  models  are  implemented  with  the  time-averaged  forms  of 
the  Navier-Stokes  equations. 

The  most  commonly  used  averaging  procedures  are: 

•  the  standard  time  averaging  procedure  for  incompressible  flow  and 

•  the  mass-averaged  approach  for  compressible  flows. 

The  time  averaging  procedure  omits  the  high  frequency  information  of  the 
turbulence,  but  the  unsteady  mean  flow  information  is  preserved.  For  the  incompressible 
case  the  randomly-  changing  flow  variables  are  replaced  with  their  averages  plus  their 
fluctuations.  Thus  the  u  velocity  component  is  represented  as  u=u  +  u'  where  u  is  the 
mean  velocity  and  u'  is  the  fluctuation  about  the  mean.  The  governing  equations  are  time 
averaged  and  the  average  of  the  fluctuation  terms  is  set  equal  to  zero. 

In  the  compressible  flow  case  the  mass-weighted  variable  of  the  Favre  averaging 

approach  is  used.  In  this  case  u  is  represented  as  u  =  u  +  u"  where  u  is 

_  pu 
u  =szzr 
P 

Here  time  average  of  the  doubly  primed  fluctuating  quantities  is  not  equal  to  zero. 


(2.8) 


After  the  substitution  is  carried  out  for  all  of  the  fluctuating  flow  variables,  the  entire 
equation  is  time  averaged.  Next,  all  the  time  averaged  terms  that  are  doubly  primed 
and  multiplied  by  density  are  defined  to  be  zero.  For  example 

pu"  =  0  (2.9) 

The  equations  of  mean  motion  resulting  from  the  time  averaging  procedure  have 
more  unknowns  than  equations.  This  constitutes  the  closure  problem  of  turbulence. 

In  order  to  close  these  equations  a  turbulence  model  must  be  used. 

The  Baldwin-Lomax  (B  -  L)  turbulence  model  [Ref.  10]  is  a  two-layer  eddy 
viscosity  model  which  simulates  the  effect  of  turbulence  in  terms  of  the  eddy  viscosity 
coefficient  pt.  The  term  p  in  the  stress  terms  is  replaced  by  p+pt  and  the  p/Pr  in  the  heat 
flux  terms  is  replaced  with  p/Pr+pt/Prt-  The  B-L  turbulence  model  bypasses  the  need  for 
finding  the  edge  of  the  boundary  layer  by  using  vorticity  instead  of  the  boundary  layer 
thickness.  This  model  is  adequate  for  flows  which  have  mild  pressure  gradients,  but  it  is 
not  very  suitable  for  highly  separated  flows.  The  basic  equations  of  the  model  follow. 

In  the  inner  layer,  the  eddy  viscosity  is  assumed  to  be  proportional  to  the  mixing 
length  squared  and  vorticity,  and  in  the  outer  layer  it  uses  an  exponentially  decaying 
formula.  The  inner  eddy  viscosity  is  computed  up  to  the  point  where  it  is  equal  to  the 
outer  eddy  viscosity  as  shown  below. 


(Pt )  inner  for 

y  —  y  crossover  I 
(ft  )  outer  for  y> 

y  crossover  J 

where  y  is  the  normal  distance  from  the  wall  and  ycrossover  taken  at  its  minimum  value 
where  it  equals  y.  The  inner  eddy  viscosity  is  given  by: 

Winner  =P^O)\ 


where 


i=*y 


1— exp(— 


\a>\= 


dco  du  2 
V  dx  dz 


7 


A+  is  an  experimentally  determined  damping  constant,  k  is  the  Von  Karman  constant. 
The  outer  eddy  viscosity  is  given  by 


III.  NUMERICAL  FLOW  SOLUTIONS 


A.  INCOMPRESSIBLE  INVISCID  FLOW 

The  basic  governing  equation  for  two-dimensional,  incompressible,  irrotational 

32<I>  92<I> 

inviscid  flow  is  the  Laplace  equation  — r-  +  — r  =  0 .  This  comes  from  the  continuity 

dx2  dz 

Eq.(2.2)  for  incompressible  flow  (p=const)  and  the  equation  for  irrotational  flow 

du  _dw 
dz  dx 

A  computer  code  developed  by  Teng  [Ref.  15]  was  used  for  the  flow  calculations. 
This  panel  code  uses  the  assumption  that  the  airfoil  can  be  approximated  by  a  number  of 
panels  each  having  a  source  strength  qj.  All  the  panels  are  assumed  to  have  the  same 
vorticity  strength  y.  These  sources  and  vortices  induce  a  velocity  at  every  point  around 
the  airfoil.  For  the  flow  solution  to  be  found  two  conditions  are  applied: 

•  The  flow  tangency  condition  which  states  that  the  normal  component  of  the 
total  velocity  at  the  midpoint  of  each  panel  due  to  ffeestream  velocity  and  the  velocities 
induced  by  the  sources  and  vortices  on  all  the  panels  is  zero. 

•  The  Kutta  condition  which  states  that  the  pressure  on  the  upper  and  lower 
surface  of  the  trailing  edge  must  be  equal. 

For  an  oscillating  airfoil  an  additional  wake  panel  is  assumed  to  be  attached  to  the 
trailing  edge.  The  airfoil  motion  is  divided  into  small  steps.  At  each  step  the  vorticity  of 
the  wake  panel  is  assumed  to  be  concentrated  into  a  single  point  vortex  which  detaches 
from  the  airfoil  with  the  local  velocity  such  that  the  Helmholtz  theorem  can  be  applied: 
the  change  in  circulation  about  the  airfoil  between  time  steps  k-1  and  k  is  equal  and 
opposite  in  direction  to  the  vorticity  released  into  the  wake  or 

=  F/t-1 
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where 


A  is  the  wake  panel  length 
yw  is  the  vorticity  strength  and 
T  is  the  circulation  about  the  airfoil 

Two  more  equations  are  needed  in  order  to  determine  the  length  and  the 
orientation  of  the  wake  panel.  Thus  two  additional  conditions  must  be  specified: 

The  wake  panel  is  oriented  in  the  direction  of  the  local  resultant  velocity  at  its  mid-point 
The  length  of  the  wake  panel  is  proportional  to  the  magnitude  of  the  local  resultant 
velocity  at  its  mid-point  and  the  time-step.  Figure  3-1  shows  the  unsteady  wake  and  its 
elements. 


A.  V  (Ik-2  “Ik-1 ) 


(Ik-3  “Ik-2) 


Figure  3-1  Incompressible  Flow  wake  model 


B.  COMPRESSIBLE  FLOW 

The  governing  equations  of  the  compressible  flow  over  airfoils  are  so  complicated 
that  no  exact  solution  can  be  obtained  analytically.  Their  solution  is  accomplished  with 
numerical  methods  which  were  applied  on  a  mass  scale  after  the  advent  of  high  speed 
computers. 
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The  total  process  of  obtaining  flow  solutions  about  problems  involving  fluid 
motion  can  be  represented  schematically  in  Figure  3-2. 


FOR  EACH  ELEMENT  OF  FLUID 


Conservation  of  mass  =>  Continuity  Equation 
Newton’s  second  law  =>  Euler  Equations 

Navier-Stokes  Equations 
Conservation  of  energy  =>  Energy  Equation 


Solve  the  equations+B.C’s+LC’s 


r 

Velocity  Distribution  : 

u(x,z,t),  w(x,z,t) 

Pressure  : 

p(X,Z,t) 

Density  : 

p(x,z,t) 

Temperature  : 

T(x,z,t) 

r 

Deduce  flow  behavior 

flow  separation 
flow  rates 
heat  transfer 
forces  on  bodies 
(skin  friction,  drag,  lift) 

Figure  3-2  Overview  of  Computational  Fluid  Dynamics 


The  governing  partial  differential  equations,  are  replaced  with  systems  of 
algebraic  equations,  so  that  a  computer  can  be  used  to  obtain  the  solution.  The  process  of 
converting  the  continuous  governing  equations  to  a  system  of  algebraic  equations  is 
known  as  discretization. 
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For  finite  differencing  schemes,  a  grid  of  discrete  points  is  distributed  throughout 
the  computational  domain  in  time  and  space.  The  algebraic  equations  link  together  values 
of  the  dependent  variables  at  adjacent  grid  points.  The  required  number  of  grid  points  for 
an  accurate  solution  typically  depends  on  three  factors;  the  dimensionality,  the  geometric 
complexity  and  the  severity  of  the  gradients  of  the  dependent  variables.  At  each  grid 
point  each  dependent  variable  and  certain  auxiliary  variables  must  be  stored.  All  of  these 
variables  must  be  stored  in  main  memory  for  the  computation  to  be  efficient. 

As  the  governing  equations  for  flows  around  an  airfoil  are  nonlinear,  the  process 
of  the  computational  solution  must  be  iterative.  The  iterative  process  is  often  equivalent 
to  advancing  the  solution  over  a  small  time  step.  The  number  of  iterations,  or  time  steps, 
might  vary  from  a  few  hundred  to  several  thousand. 

The  discretization  process  always  introduces  an  error.  As  long  as  the  discrete 
equations  are  faithful  representations  of  the  governing  equations,  this  error  can  be 
reduced  by  refining  the  grid.  If  the  numerical  algorithm  that  performs  the  iteration  is 
stable,  then  the  computational  solution  can  be  made  arbitrarily  close  to  the  true  solution 
of  the  governing  equations  by  refining  the  grid. 

1 .  Grid  Generation  by  Algebraic  Mapping 

In  order  to  use  an  unweighted  differencing  scheme  the  flow  equations  must  be 
transformed  to  a  generalized  coordinate  system  using  the  following  transformations  for 
the  case  of  unsteady  two-dimensional  flow: 

4  =  £(x,z,t) 

4  =C(x,z,t)  (3.1) 

t=z(x,z,t) 

The  above  equations  can  be  used  to  transform  the  governing  equations  from  the 
physical  (x,  z,t)  to  the  computational  domain  (£,,^,x)  with  the  Jacobian  transformation. 

For  more  details  refer  to  [Ref. 2]. 
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Grid  generation  has  to  do  with  the  establishment  of  the  correspondence  between 
points  (x,  z)  in  the  physical  domain  and  points  (£,Q  in  the  computational  domain. 

Algebraic  mapping  techniques  interpolate  the  boundary  data  in  one  or  more 
dimensions  in  order  to  generate  the  interior  grid.  The  generated  grid  should  be  well- 
conditioned,  i.e.,  smoothly  varying,  close  to  orthogonal  and  with  local  grid  aspect  ratios 
close  to  unity.  For  fluid  flow  problems,  the  solution  is  often  changing  rapidly  close  to  a 
particular  surface.  It  is  important  to  construct  a  grid  that  is  orthogonal,  or  near- 
orthogonal,  adjacent  to  such  a  surface. 

Stretching  functions  on  the  boundaries  are  used  in  order  to  define  the  grid  points 
in  the  interior.  Consequently,  the  two-boundary  and  multisurface  techniques  are  still  able 
to  get  smoothly  varying  near-orthogonal  grids  with  only  one-dimensional  explicit 
interpolation. 

The  distribution  of  points  along  the  boundary  of  the  domain  is  handled  effectively 
by  defining  normalized,  one-dimensional,  stretching  functions  along  boundary  segments, 
typically  corresponding  to  each  side  of  the  computational  rectangle  in  the  (£77)  plane. 
Boundary  stretching  functions  are  applicable  whether  the  interior  grid  is  generated  by 
solving  a  partial  differential  equation  or  by  an  algebraic  mapping. 

2.  Numerical  Implementation  of  Algebraic  Mapping 

Computer  programs  have  been  developed  capable  of  generating  grids  between 
two  bounding  curves  based  on  algebraic  mapping  techniques. 

For  the  domain  shown  in  (Figure  3-3),  the  bounding  surfaces  consist  of  a 
symmetric  slender  body  extended  downstream,  ABC,  and  a  farfield  boundary,  FED. 
Between  these  two  boundaries  half  of  a  C-grid  is  to  be  generated.  By  symmetry  the 
complete  C-grid  can  be  obtained  by  reflection  about  the  x-axis. 

As  mentioned  before,  grid  generation  is  split  into  two  parts: 

First,  grid  point  locations  on  all  boundaries  are  determined  and  a  1-D  stretching 
function  is  used  to  control  the  distribution  on  the  boundaries. 
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Subsequently,  the  interior  grid  is  generated  by  the  multisurface  technique.  Two 
intermediate  surfaces,  Z2,  and  Z3,  are  introduced,  one  each  adjacent  to  the  bounding 


surfaces  ABC  and  FED.  The  parametric  (r)  correspondence  of  surfaces  Z2,  and  Z3  to  their 
neighbouring  bounding  surface  is  adjusted  so  that  grid  lines  intersect  the  bounding 
surfaces  orthogonally.  The  mechanism  of  choosing  x(r),  y(r)  on  surfaces  Z2,  and  Z3 
requires  an  orthogonal  projection,  conceptually  similar  to  the  near-orthogonal  grid 
constmction. 


3.  CFD  Techniques 

Computational  techniques  are  used  to  obtain  an  approximate  solution  of  the 
governing  equations  and  boundary  conditions.  For  example,  for  two-dimensional 
unsteady  incompressible  flow,  velocity  and  pressure  solutions,  u(x,  z,  t),  w(x,  z,  t)  and 
p(x,  z,  t),  would  be  computed.  The  computational  solution  is  obtained  in  two  steps  that 
are  shown  in  Figure  3-4. 

In  the  first  step,  the  partial  differential  equations  and  boundary  and  initial 
conditions  are  converted  into  a  discrete  system  of  algebraic  equations.  This  step,  as 
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mentioned  before,  is  called  discretization.  The  second  step  comprises  the  solution  of  the 
system  of  algebraic  equations. 

The  fact  that  the  differentiated  terms  in  the  governing  partial  differential 
equations  are  replaced  by  algebraic  expressions  connecting  nodal  values  on  a  finite  grid 
introduces  an  error.  In  order  to  minimize  this  error,  the  most  appropriate  algebraic 
expressions  should  be  chosen.  Equally  important  as  the  error  in  representing  the 
differentiated  terms  in  the  governing  equation  is  the  error  in  the  solution. 


Figure  3-4  Overview  of  the  computational  solution  procedure 


a.  Discretisation  process  example 

To  convert  the  governing  partial  differential  equation(s)  into  a  system  of 
algebraic  equations  (or  ordinary  differential  equations),  a  number  of  choices  are 
available.  The  most  common  are  the  finite  difference,  finite  element,  finite  volume  and 
spectral  methods. 

The  way  the  discretisation  is  performed  also  depends  on  whether  time 
derivatives  (in  time  dependent  problems)  or  equations  containing  only  spatial  derivatives 
are  being  considered.  In  practice,  time  derivatives  are  discretized  almost  exclusively 
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using  the  finite  difference  method.  Spatial  derivatives  are  discretized  by  either  the  finite 
difference,  finite  element,  finite  volume  or  spectral  method. 

The  discretization  process  can  be  illustrated  by  considering  the  scalar 

equation 

du  du 

— +  «— =  °  (3.2) 

ot  ox 

The  most  direct  means  of  discretization  is  provided  by  replacing  the 
derivatives  by  equivalent  finite  difference  expressions.  It  is  based  on  forward  difference 
in  time  and  a  central  difference  in  space  and  gives,  after  multiplying  with  At: 

«r1-<,+frte.-«w)= 0  (3-3) 

2Ax 

This  equation  is  referred  to  as  explicit  in  time.  That  is  if  we  know  u"  at  all 
grid  points  in  space  at  time  level  n,  then  we  have  a  series  of  explicit  calculations  to 
determine  un+1  at  all  grid  points.  Further,  all  values  of  Uin+1  are  obtained  independently  of 
each  other;  they  depend  only  on  u",  not  on  the  other  Uin+1  terms. 

The  process  of  discretizing  Eq.(3.2)  to  give  Eq.(3.3)  implies  that  the 
problem  of  finding  the  exact  (continuous)  solution  u(x,  t)  has  been  replaced  with  the 
problem  of  finding  discrete  values  Uj“,  i.e.  the  approximate  solution  at  the  (i,  n)th  node.  In 
turn,  two  related  errors  arise,  the  truncation  error  and  the  solution  error. 

The  precise  value  of  the  approximate  solution  between  the  nodal  (grid) 
points  is  not  obvious.  Intuitively,  the  solution  would  be  expected  to  vary  smoothly 
between  the  nodal  points.  In  principle,  the  solution  at  some  point  (xr,  ur )  that  does  not 
coincide  with  a  node  can  be  obtained  by  interpolating  the  surrounding  nodal  point 
solution. 

To  provide  the  complete  numerical  solution  at  time  level  (n  +  1),  Eq.  (3.3) 
must  be  applied  for  all  the  nodes  i=2,...I-l,  assuming  that  Dirichlet  boundary  conditions 
provide  the  values  uin+1  and  uin+1. 

The  discretization  process  invariably  introduces  an  error  unless  the 
underlying  exact  solution  has  a  very  elementary  analytic  form.  In  general,  the  error  for  a 
finite  difference  representation  of  a  derivative  can  be  obtained  by  making  a  Taylor  series 
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expansion  about  the  node  at  which  the  derivative  is  being  evaluated.  The  evaluation  of 
the  leading  term  in  the  remainder  provides  a  close  approximation  to  the  error  if  the  grid 
size  is  small.  However,  the  complete  evaluation  of  the  terms  in  the  Taylor  series  relies  on 
the  exact  solution  being  known. 

b.  Convergence 

A  solution  of  the  algebraic  equations  that  approximate  a  given  partial 
differential  equation  is  said  to  be  convergent  if  the  approximate  solution  approaches  the 
exact  solution  of  the  partial  differential  equation  for  each  value  of  the  independent 
variable  as  the  grid  spacing  tends  to  zero.  Thus  we  require 

m"  — » u(xi,tn )  as  Ax, At  — > 0 
where  u  is  the  exact  solution 

The  difference  between  the  exact  solution  of  the  partial  differential 
equation  and  the  exact  solution  of  the  system  of  algebraic  equations  is  the  solution  error, 
denoted  by  e;  that  is 

e"  =M(x(.,r„)-M,n 

The  exact  solution  of  the  system  of  algebraic  equations  is  the  approximate  solution  of  the 
governing  partial  differential  equation.  The  exact  solution  of  the  system  of  algebraic 
equations  is  obtained  when  no  numerical  errors  of  any  sort,  such  as  those  due  to  round¬ 
off,  are  introduced  during  the  computation.  The  magnitude  of  the  error,  at  the  (bn)*  node 
typically  depends  on  the  size  of  the  grid  spacings,  Ax  and  At,  and  on  the  values  of  the 
higher-order  derivatives  at  that  node,  omitted  from  the  finite  difference  approximations  to 
the  derivatives  in  the  given  differential  equation. 

4.  Navier  -  Stokes  Solution 

If  we  apply  the  generalized  transformation  to  the  compressible  Navier-Stokes 
equations  written  in  vector  form  Eq.(2.7),  the  following  transformed  equation 
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is  obtained 


bq  a#  3g=j_  as 
dt+d£  +  d£  Red£ 

where  Q  is  the  conservative  variables  vector: 

P 

pu 

pw 
e 

J  is  defined  in  Eq.(3.9) 

F  ,G  are  the  in  viscid  flux  vectors: 

pU  1  pW 

puU  +£xp  ^  =  J_  /wW  +  £xp 

pwU  +  tjfP  J  pwU  +  £zp 

(e  +  p)u  -  £t  p\  |_0  +  p)W  -  £  p_ 

and  S  is  the  thin  layer  approximation  of  the  viscous  fluxes  in  the  £  direction: 

0 

pmxU£  +(///3 )m2£x 
pmxW£  +  (p/3)m2£z 
pmxm-i  +  (p/3  )m2m4 

where 

ro3  =(m2  +w2)/2+  ^ —  } 

3  Pr(r-i) 

«4=£,“  +  £zW 


(3.4) 


The  U  and  W  are  the  contravariant  velocity  components,  tangential  to  the  constant  £,  and 
^  surfaces,  respectively,  and  they  are  given  by: 
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u=u4,  +  w£t  +4, 

W=u{,+w(,+(, 


Non-dimensionalization  of  the  governing  Eq.(3.4)  is  performed  using  the 
following  reference  quantities:  density  p~,  length  c,  time  c/a^,  and  energy  p„a2„ 

After  non-dimensionalization  of  the  governing  equations,  Euler  solutions  are 
obtained  setting  the  viscous  terms  of  Eq.(3.4)  zero  and  applying  the  flow  tangency 
condition  at  the  airfoil  surface.  For  Navier-Stokes  solutions  the  no-slip  condition  is 
applied  at  the  surface. 

The  numerical  integration  of  the  Eq  (3.4)  is  performed  using  an  implicit  scheme 
[Ref.  1 1]  given  by 


|A*,+A;A-J 


(3.5) 


1.0 Uc  Qi,k  ^i-V2,k  )+  (g,>+1/2  ^a-1/2  )  ^  (^++1/2  —  Si,k-U2  )J 


In  Eq.(3.5)  V,Aand  8  are  the.  forward,  backward  and  central  difference  operators 
respectively,  whereas  the  hi=A'i/A<!; ,  h(  =At  I  At,  ,  A±  =dF/dQand 

B±  =dG/dQ  are  the  flux  Jacobian  matrices.  The  superscript  (.)”  refers  to  the  time  step 
and  the  superscript  (.)p  refers  to  Newton  subiterations  within  each  time  step. 

Time  accuracy  of  the  implicit  numerical  solution  is  obtained  by  performing 
Newton  iteration  to  convergence  for  each  time  step.  Linearization  and  factorization  errors 
are  minimized  because  the  left  hand  side  of  Eq.(3.5)  can  be  driven  to  zero  at  each  time 
step.  Typically  two  to  three  subiterations  are  sufficient  to  drop  the  residuals  two  orders  of 
magnitude  during  the  Newton  iteration  process. 
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IV.  AEROELASTIC  ANALYSIS 

A.  FLUTTER  -  GENERAL 

Flutter  can  be  defined  as  the  dynamic  instability  of  an  elastic  body  in  an  airstream. 
The  aircraft  aerodynamic  surfaces,  such  as  wings,  tails,  and  control  surfaces  (i.e.  bodies 
subjected  to  large  lateral  aerodynamic  loads  of  the  lift  type)  are  the  bodies  that  are  most 
likely  to  experience  flutter.  The  only  air  forces  necessary  to  produce  it  are  those  due  to 
deflections  of  the  elastic  structure  from  the  undeformed  state. 

The  flutter  speed  Up  and  frequency  (0F  are  defined,  respectively,  as  the  lowest 
airspeed  and  the  corresponding  circular  frequency  at  which  a  given  structure  flying  at 
given  atmospheric  density  and  temperature  will  exhibit  sustained,  simple  harmonic 
oscillations.  All  small  motions  are  stable  at  speeds  below  Uf,  whereas  divergent 
oscillations  can  occur  in  a  range  of  speeds  (or  even  at  all  speeds)  above  Uf.  In  other 
words,  flight  at  Up  represents  a  borderline  condition  or  neutral  stability  boundary  above 
which  flutter  will  happen. 

.  The  free  vibration  of  a  linear  structure  in  vacuum  is  a  real,  or  single,  eigenvalue 
problem.  In  contrast,  the  theoretical  flutter  analysis  leads  to  a  complex,  or  double, 
eigenvalue  problem,  where  two  characteristic  numbers  determine  the  critical  flutter  speed 
and  frequency.  This  is  because  the  simple  harmonic  motion  assumption  is  made;  all 
dependent  variables  are  proportional  to  eI<Bt,  and  the  problem  is  to  find  all  the 
combinations  of  U  and  co  for  which  flutter  actually  occurs. 

There  are  three  main  ways  in  order  to  study  and  predict  flutter: 

•  theoretical  flutter  computation 

•  wind-tunnel  experiments  on  scaled  dynamic  models, 

•  flight  testing  of  full-scale  aircraft. 
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The  decision  as  to  which  of  these  is  most  economical  in  a  given  case  is  very  difficult  and 
depends  on  a  multitude  of  factors  such  as: 

•  the  purpose  and  the  phase  of  the  study 

•  anticipated  margin  of  safety  from  flutter, 

•  the  Mach  number  range, 

•  the  number  of  different  mass  and  structural  configurations  to  be  analyzed. 
Current  regulations  require  that  flight  flutter  testing  must  be  preceded  by  flutter  analysis 
and  wind  tunnel  testing. 


B.  THE  PHENOMENON  OF  FLUTTER 

The  phenomenon  of  flutter  can  be  described  with  a  simple  experiment:  consider  a 
cantilever  wing,  mounted  in  a  wind  tunnel  at  a  small  angle  of  attack  and  with  the  root 
rigidly  built  in.  Then  when  there  is  no  flow  in  the  tunnel  we  deliberately  deflect  and  then 
release  the  wing.  We  observe  that  the  oscillation  decays  gradually.  Then  we  increase  the 
speed  of  flow  in  the  wind  tunnel  and  repeat  the  experiment;  the  oscillation  of  the 
disturbed  airfoil  dies  out.  With  further  increase  of  the  flow  speed,  however,  a  point  is 
reached  at  which  the  oscillation  maintains  itself  with  steady  amplitude.  This  is  the  critical 
flutter  speed.  If  the  wind  tunnel  speed  is  increased,  flutter  will  occur  with  increasing 
amplitude.  It  is  obvious  that  during  the  aircraft  flight  at  speeds  above  the  critical,  an 
initial  disturbance  (coming  from  a  gust  for  example)  will  cause  an  oscillation  of  great 
violence.  In  such  circumstances  the  airfoil  suffers  from  oscillatory  instability  and  is  said 
to  flutter. 

The  above  experiment  on  wing  flutter  shows  that  the  oscillation  is  self-excited; 
i.e.,  no  further  external  disturbance  is  required.  The  motion  can  maintain  itself  or  grow 
for  a  range  of  wind  speeds,  which  depends  on  the  design  of  the  wing  and  the  conditions 
of  the  test. 
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An  aerodynamic  surface,  as  an  elastic  body,  has  infinitely  many  degrees  of 
freedom  (d.o.f).  But  often  it  is  sufficient  to  consider  only  three  variables,  namely  the 
flexure  (plunge),  the  torsion  (pitch),  and  the  control-surface  rotation.  A  flutter  mode 
consisting  of  all  three  elements  is  called  a  ternary  flutter.  In  special  cases,  however,  two 
of  the  variables  predominate,  and  the  corresponding  flutter  modes  are  called  binary  flutter 
modes.  In  fact,  many  airplanes  can  be  replaced  by  a  system  of  simple  beams,  so  that  the 
elastic  deformation  can  be  described  by  the  deflection  and  torsion  of  the  beams,  in 
addition  to  the  rotation  of  control  surfaces  about  their  hinge  lines. 

For  the  oscillatory  motion  of  a  wing  only  the  first  two  components  of  flutter  are 
often  considered:  the  flexural  (plunge)  and  the  torsional  (pitch).  In  general,  the  coupling 
of  those  two  degrees  of  freedom  is  an  essential  feature  for  wing  flutter.  The  oscillation 
that  occurs  at  the  critical  speed  is  harmonic.  Experiments  on  cantilever  wings  show  that 
the  flexural  movements  at  all  points  across  the  span  are  approximately  in  phase  with  one 
another,  and  likewise  the  torsional  movements  are  all  approximately  in  phase,  but  the 
flexure  is  considerably  out  of  phase  from  the  torsional  movement.  As  will  be  seen  later,  it 
is  this  phase  difference  that  is  responsible  for  the  occurrence  of  flutter. 

A  rigid  airfoil  in  a  low-speed  attached  flow  which  is  constrained  in  torsion  and 
can  only  plunge  (flexural  d.o.f)  does  not  flutter.  In  contrast,  a  rigid  airfoil  with  only  the 
torsional  d.o.f  can  flutter  for  some  special  mass  distributions  and  elastic-axis  locations.  In 
this  text  we  will  restrict  the  term  “flutter”  to  the  oscillatory  instability  in  a  flow  which 
does  not  experience  flow  separation  or  strong  shocks. 

For  a  control  surface,  such  as  a  flap  or  an  aileron,  its  flexural  and  torsional  elastic 
deformation  are  not  so  important  as  its  freedom  to  turn  about  the  hinge  line. 

Consequently,  the  deflection  of  a  control  surface  can  simply  be  described  by  the  angle  of 
rotation  about  its  hinge  line. 

The  degrees  of  freedom  of  flutter  mentioned  before,  together  with  the  freedom  of 
the  airplane  to  move  as  a  rigid  body,  offer  a  large  number  of  possible  combinations  of 
binary,  ternary,  and  higher  modes  of  flutter.  Since  it  is  not  clear  which  of  these  modes 
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correspond  to  the  actual  critical  speeds,  it  is  necessary  either  to  resort  to  experiments  or, 
in  a  theoretical  approach,  to  analyze  all  cases.  This  is  why  a  successful  flutter  analysis 
depends  so  much  on  the  analyst’s  experience.  He  must  be  able  to  choose,  among  all 
possible  modes,  those  that  are  likely  to  be  critical  for  a  given  structure. 

As  one  can  easily  understand,  flutter  analysis  is  an  extensive  subject;  for  this 
reason  it  is  preferable  to  divide  its  presentation  in  sections.  Some  basic  definitions  are 
given  in  Section  C.  The  role  of  the  elastic  stiffness  in  flutter  prevention  are  explained  on 
an  empirical  basis  in  Section  D.  The  origin  of  flutter  from  the  aerodynamic  point  of  view 
is  then  considered  in  Section  E.  It  will  be  shown  that  flutter  occurs  because  the  speed  of 
flow  affects  the  amplitude  ratios  and  phase  shifts  between  motion  in  various  degrees  of 
freedom  in  such  a  way  that  energy  can  be  absorbed  by  the  airfoil  from  the  airstream 
passing  by.  The  physical  explanation  and  the  equations  of  motion  for  the  one  d.o.f  flutter 
will  be  discussed  in  Section  F,  whereas  the  two  d.o.f  flutter  will  be  explained  in  Section 
G.  A  practical  way  to  find  the  flutter  condition  with  the  Panel  Code  (UPOT)  will  be 
shown  in  Section  H. 


C.  NONDIMENSIONAL  PARAMETERS 

In  flutter  analysis  any  non-dimensional  quantity  relating  to  the  motion  can  be 
expressed  as  a  function  of  two  parameters: 

p/a  and  K/ai3U2 

where 

1  Typical  linear  dimension 

U  Air'  speed 

p  Air  density 

a  Typical  density  of  structural  material 

K  Typical  torsional  stiffness  constant  (ft-lb  per  rad) 
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If  we  consider  the  energy  dissipation  of  the  structure,  the  viscosity  and  the 
compressibility  of  the  fluid  in  flutter  analysis  then  three  more  non-dimensional  variables 
must  be  added: 

g  Material  damping  coefficient 

Re  Reynolds  number 

M  Mach  number 

For  two  systems  to  be  dynamically  similar  they  should  have  equal  variables  g,  Re, 
and  M  in  addition  to  the  above  parameters.  In  general,  g  is  important  in  control-surface 
flutter.  Re  is  important  in  stall  flutter,  and  M  is  important  in  high-speed  flight;  otherwise 
their  effects  are  small. 

The  non-dimensional  form  of  the  frequency  of  oscillation  00  (radians  per  second, 
with  dimension  T'1)  is: 


which  is  called  the  reduced  frequency. 

Two  similar  systems  having  the  same  values  of  p/<7and  K/al3U2  flutter  at  the 
same  reduced  frequency.  Since  all  derived  concepts  relating  to  the  motion  can  be 
expressed  in  functional  relations  as  above,  it  is  clear  that  the  equality  of  the  values  of  the 
parameters  p/crand  K/al3U2  is  sufficient  to  guarantee  dynamic  similarity  of  the  two 
systems. 

The  reduced  frequency  characterizes  the  variation  of  the  flow  with  time.  Its 
inverse,  U I  col  is  called  the  reduced  speed.  The  physical  meaning  of  the  reduced 
frequency  is  given  as  follows:  Consider  that  a  disturbance  occurs  at  a  point  on  a  body  and 
oscillates  together  with  the  body.  The  fluid  element  influenced  by  the  disturbance  moves 
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downstream  with  a  mean  velocity  U.  Let  the  frequency  of  oscillation  of  the  body  and  the 
disturbance  be  co.  Then  the  spacing,  or  "wave  length"  of  the  disturbance,  is  X=  T-U  or 
X=2nU/(0.  So  the  reduced  frequency  for  the  airfoil  of  chord  c  is: 


,  coc  .  c 
k=  —  =  In— 
U  X 


(4.2) 


Therefore,  the  reduced  frequency  represents  a  ratio  of  the  characteristic  length  of  the  body 
(chord)  to  the  wave  length  of  the  disturbance.  In  other  words,  the  reduced  frequency 
characterizes  the  way  a  disturbance  is  felt  at  other  points  of  the  body.  Since  every  point 
of  an  oscillating  body  disturbs  the  flow,  one  may  say  that  the  reduced  frequency 
characterizes  the  mutual  influence  between  the  motion  at  various  points  of  the  body. 

The  reduced  natural  pitching  frequency  is  defined  as 

CQ„c 


k„  =- 


U 


(4.3) 


where 


(4.4) 


is  the  uncoupled  natural  torsional  frequency  of  the  system,  K0is  the  spring  constant  for 
pitching  and  Ia  is  the  moment  of  inertia  about  the  pivot  point  per  unit  span.  From  Eq.(4.3) 
and  (4.4)  we  have: 


(4.5) 


The  reduced  flutter  velocity  is  defined: 


U 

0)ac 


(4.6) 
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D.  STIFFNESS  CRITERIA 


As  mentioned  before,  wing-flutter  occurs  when  the  reduced  frequency  is  lower 
than  the  following  critical  value: 


k 


cr 


U 


where  U  is  the  mean  speed  of  flow,  coa  is  the  undamped  natural  pitching  frequency  of  the 
wing,  and  c  is  the  chord  length  of  the  vibrating  portion  of  the  wing. 

For  safety  against  flutter,  the  reduced  frequency  should  be  higher  than  kcr,  or  the 
design  speed  of  the  airplane  should  be  lower  than 


v„-¥- 

k„ 


(4.7) 


The  frequency  coa,  which  may  be  determined  by  ground  vibration  experiments  or 
computed  by  a  theoretical  analysis,  increases  with  increasing  stiffness  of  the  structure. 
Therefore,  the  critical  speed  can  be  raised  by  increasing  the  wing  stiffness.  It  is  obvious 
that  when  the  lower  limit  of  Ucr  is  defined  (e.g.,  the  maximum  speed  of  flight)  we  can 
determine  the  minimum  value  of  (0,  and  hence  the  minimum  value  of  the  torsional 
rigidity. 


E.  VERTICAL  TRANSLATORY  OSCILLATION 

The  fundamental  role  played  by  the  aerodynamic  forces  in  inducing  flutter  is  that 
they  are  the  means  through  which  kinetic  energy  is  transferred  from  the  airstream  to  the 
airfoil  as  elastic  and  kinetic  energy.  Hence,  the  possibility  of  flutter  can  be  discussed  by 
considering  the  energy  relation.  From  the  analysis  we  can  find  the  amount  of  energy 
exchanged  between  the  airstream  and  the  airfoil;  if  the  oscillating  body  gains  energy  from 
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the  airstream  in  completing  a  cycle  then  the  oscillation  is  aerodynamically  unstable;  if  it 
gives  energy  to  the  airstream  the  oscillation  is  going  to  be  aerodynamically  stable. 

In  order  to  determine  the  energy  transfer  consider  an  airfoil  performing  a  vertical 
translatory  oscillation  with  a  constant  amplitude  h0.  Then  the  vertical  displacement  can 
be  described  by 

h  =  h0eim  (4.8) 

where  h  is  defined  to  be  as  positive  downward.  The  speed  of  downward  motion  is 
therefore 

h  =  ia>h0eim  '  (4.9) 

where  the  prime  indicates  a  differentiation  with  respect  to  time.  If  h  is  a  constant,  the 
downward  motion  will  induce  a  lift  force  L0  on  the  airfoil: 


h=\pu2s 


dCL  h 

HaU 


(4.10) 


where  £.~ _ is  the  dynamic  pressure  and  S  is  the  wing  area.  This  value  L0  is  the  quasi¬ 

steady  lift  and  is  defined  as  positive  upward  in  the  usual  sense.  As  the  airfoil  performs  a 
translatory  oscillation,  the  true  instantaneous  lift  acting  on  the  airfoil  differs  from  L0  both 
in  magnitude  and  in  phase  and  its  value  can  be  stated  in  the  form 

L  =  L0rei¥  (4.11) 

In  this  expression  r  represents  the  ratio  of  the  absolute  value  of  the  instantaneous  lift  to 
that  of  the  quasi-steady  lift,  and  \\r  the  phase  angle  by  which  the  actual  lift  leads  the  quasi¬ 
steady  value.  The  quantities  r  and  \\f  depend  on  the  reduced  frequency  k,  the  Mach 
number  M,  and  the  Reynolds  number  Re.  For  a  nonviscous  incompressible  fluid,  r  and  \y 
are  functions  of  k  alone.  The  ratio  L/L0  =  reiv  can  be  plotted  as  vector  with  length  r  and 
angle  \\f.  Using  incompressible  oscillatory  thin  airfoil  theory  von  Karman  and  Sears 
derived  the  vector  diagram  of  Figure  4-1. 
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Figure  4-1  Vector  diagram  of  lift  in  vertical  translation  oscillation.  (From  Y.C.Fung  Ref.  1) 


When  the  airfoil  moves  through  a  distance  dh,  the  differential  work  done  by  the 

lift  is 


dW  =  -Ldh  =  -Lhdt  (4.12) 

and  therefore  the  total  work  done  by  the  air  on  the  airfoil  during  one  cycle  of  oscillation, 
as  given  by  Y.C  Fung,  [Ref.  1]  is: 


2tcI(o  1  2iz!(o 

W  =  -  (  Re[L]-Re[h]  dt  =  -—pUS—L(coh0)2r  f  sin(cot  +  ys)  ■  sin(ax)dt  = 
«  2  da  i 


K  dC,  o 

= — pUS — —coh0  rcosy/ 
2  da 


(4.13) 


Hence,  the  gain  of  energy  Wby  the  airfoil  from  the  airstream  is  proportional  to 
-cos\|/,  where  \jr  is  the  phase  angle  by  which  the  actual  lift  leads  the  quasi-steady  value.  If 
-Jt/2  <\|t  <  Tjfl,  then  W  is  negative;  i.e.,  the  oscillating  airfoil  will  give  energy  to  the 
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airstream  and  the  oscillation  is  therefore  stable.  If  we  refer  back  to  Figure  4-1,  it  is  seen 
that  the  condition  -Jt/2  <\j/  <  jt/2  is  satisfied.  Hence,  in  a  nonviscous  incompressible 
fluid,  the  vertical  translation  oscillation  is  aerodynamically  stable. 

This  example  shows  the  importance  of  the  phase  angle  between  the  aerodynamic 
force  and  the  oscillatory  motion.  Based  on  this  phase  angle  we  can  conclude  that  purely 
translational  flutter  in  one  degree  of  freedom  is  impossible. 

However  as  shown  in  the  next  Section,  this  conclusion  does  not  hold  for  single- 
degree-of-freedom  torsional  motion. 

F.  ONE  DEGREE  OF  FREEDOM  FLUTTER  (TORSIONAL) 

In  this  section  one  degree  of  freedom  torsional  flutter  will  be  examined.  In 
Subsection  1  the  negative  damping  at  low  reduced  frequency  k  will  be  explained 
qualitatively.  In  Subsection  2  the  time  and  frequency  domain  analysis  will  be  presented. 

1 .  Physical  Explanation 

•  Consider  an  airfoil  in  an  incompressible  airflow  of  speed  U,  subjected  to  a  slow 
(low  kj,  nearly  sinusoidal  torsional  motion.  Assume  that  the  rotational  axis  is  at  the 
leading  edge  (L.E.),  but  similar  results  are  obtained  if  it  is  anywhere  in  front  of  the  quarter 
chord  point. 

Examine  now  step  by  step  how  flutter  may  occur  for  this  one  degree  of  freedom 
oscillating  airfoil  with  the  help  of  Figure  ^  given  by  Bisplinghoff  et  al  [Ref. 4].  At  stage 
(a)  when  the  airfoil  is  at  the  highest  angle  oi  attack  all  the  vortices  shed  during  the 
preceding  last  half  cycle  have  counter-clockwise  vorticity.  Between  (a)  and  (c)  the  airfoil 
pitches  down  resulting  in  vortices  of  clockwise  vorticity.  As  k  is  low,  the  counter¬ 
clockwise  vortices  are  far  away  from  the  trailing  edge  (T.E.)  and  do  not  contribute 
significantly  in  inducing  the  upwash.  Therefore,  as  the  airfoil  pitches  through  the  zero 
AO  A  position,  stage  (b),  a  lift  L2  is  induced  which  is  in  the  same  direction  as  the  motion 
and  therefore  net  positive  work  is  done  on  the  airfoil. 
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When  the  airfoil  increases  the  angle  of  attack  from  steps  (c)  to  (a)  the  counter¬ 
clockwise  vortices  dominate  inducing  a  wake  upwash  and  therefore  a  lift  L2  and  moment 
which  is  in  the  same  sense  as  the  angular  velocity  and  again  does  positive  work  on  the 
wing. 

This  transfer  of  energy  from  the  airstream  to  the  airfoil  is  the  origin  of  the 
negative  damping.  At  higher  values  of  k  the  wave  length  of  the  vertical  wake  becomes 
shorter  so  that  vortices  from  previous  half-cycles  can  contribute.  The  net  effect  is  an 
eventual  change  of  the  aerodynamic  damping  from  negative  to  positive. 
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The  phenomenon  of  upwash  and  lift  induced  by  the  wake  vortices  can  be 
visualized  by  the  incompressible  panel  code  UPOT.  The  airfoil  wake  vortices  are  shown 
in  Figure  4-3  as  squares.  The  solid  ones  are  the  clockwise  vortices  whereas  the  dotted 
ones  are  the  counterclockwise  vortices.  The  airfoil  is  shown  in  the  same  four  positions  as 
in  Figure  4-2.  It  is  clearly  seen  that  the  rotation  of  the  vortices  sketched  by  Bisplinghoff  et 
al  is  reproduced  by  the  panel  code. 


Figure  4-3  Incompressible  flow  wake  visualization  (from  UPOT) 
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The  essence  of  the  foregoing  discussion  is  the  inclusion  of  unsteady  effects.  If 
only  quasi-steady  aerodynamics  had  been  used  (without  vortex  shedding),  the  possibility 
of  torsional  flutter  would  have  been  overlooked.  For  a  more  detailed  analysis  refer  to 
Bisplinghoff  et  al  [Ref.  4  pages  262-263], 

2.  Flutter  Analysis 

Flutter  analysis  can  be  done  with  two  methods.  The  first  is  the  time-domain 
analysis  in  which  the  airfoil  motion  is  found  as  a  function  of  time.  It  is  based  on 
discretisation  of  the  ordinary  differential  equation  of  motion.  The  second  method  is  the 
classical  theoretical  approach  based  on  frequency-domain  analysis  with  which  the  critical 
flutter  speed  and  frequency  are  calculated.  For  both  methods  we  consider  a  simple  system 
of  a  two-dimensional,  rigid  airfoil  of  unit  span,  which  is  hinged  at  a  specific  point  (pivot) 
but  elastically  restrained  from  rotating  about  this  axis  by  a  torsion  spring  with  constant 
Ka.  The  airfoil,  which  has  a  moment  of  inertia  about  the  pivot  point  I0,  is  placed  in  an 
airstream  so  that  the  unstrained  position  of  the  spring  corresponds  to  zero  angle  of  attack 
a  (Figure  4-4).  Then  we  let  the  airfoil  perform  a  torsional  oscillation. 
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a  -  s 


(4.18) 


(4.19) 


This  is  the  system  of  equations  used  for  the  incompressible  flutter  analysis 
where  cm  is  computed  at  each  time  step  using  the  panel  code  UPOT. 

For  the  compressible  flutter  analysis  the  time  can  be  non-dimensionalized 
using  the  free  stream  speed  of  sound,  a„ ,  instead  of  the  freestream  velocity  U „ .  Thus  the 
dimensionless  time  can  be  defined  as 


to 

t  =  —  (4.20) 

c 

Using  this  non-dimensional  cm,  ia  and  x,  the  equation  of  motion  (4.19)  becomes: 

s'  =  -kla.  +  —Mlcm 


The  reduced  natural  frequency  kafor  this  case  is  based  on  the  freestream 
speed  of  sound,  whereas  in  the  incompressible  case  it  is  based  on  the  freestream  velocity. 
In  order  to  have  compatibility  between  the  values  of  ka  found  for  incompressible  and 
compressible  flow  we  have  to  convert  the  latter  into  the  former  values.  The  procedure 
will  be  shown  in  Section  H. 

Integration  of  Eq.(4.18)  and  (4.19)  can  be  done  with  a  first-order  Euler 
scheme,  a  second  order  modified  Euler  scheme  or  a  fourth  order  Runge-Kutta  scheme. 

The  second-order  modified  Euler  scheme  is  derived  by  applying  the 
trapezoidal  rule  to  integrate  an  equation  y  =  f(x,t )  as  follows: 

y»+i  =  y,  +h/2  [f(yn+1,tn+1)  +  f(yn,tnj] 

and  using  an  initial  value  for  y(0). 

The  fourth-order  Runge-Kutta  scheme  is  derived  by  applying  the  following 

equation: 

yn+i  =  yn  +(l/6)fc  +  2k1+2ki+kA] 

with 
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*1  =¥(yn’tn )  ¥  =W(yn  +kl/2,tn+h/2) 

k3=hf(yn  +  k2/2,t„+h/2)  k4=hf(yn+k3,tn+h) 

Panel  codes,  like  UPOT,  typically  use  large  time  steps  and  therefore  use  of 
the  fourth-order  Runge-Kutta  scheme  is  advisable.  In  contrast,  the  stability  limits  of 
Euler/N-S  codes  typically  require  small  time  steps  and  therefore  the  first-order  Euler 
scheme  gives  satisfactory  accuracy. 


b.  Frequency  domain  analysis 


In  order  to  perform  the  classical  flutter  analysis  in  the  frequency  domain 
we  assume  that  the  airfoil  executes  a  torsional  oscillation  with  a  constant  amplitude  cxo- 
The  angle  of  attack  a  at  any  time  can  be  described  by 

a  =  a0e,m  (4.21) 

and  after  substituting  (4.16)  into  (4.14),  and  dividing  through  by  npb4  or  a^e10*  we  get 

\2  "~ 


71  pb 


i_r^o 

j 


+  C„  =0 


(4.22) 


where  coo  is  the  undamped  natural  frequency  of  torsional  vibration  and  cm  is  the 
dimensionless  aerodynamic  coefficient  which  is  a  complex  number. 

Hence  (4. 18)  can  be  split  into  real  and  imaginary  parts 

\2 

i  I  / 


71  pu 

Mpm}=o 


0) 


(4.23) 

(4.24) 


Flutter  occurs  only  at  those  values  of  the  reduced  frequency  that  make  the 
out-of-phase  component  of  the  aerodynamic  moment  Eq.(4.24)  vanish,  provided  that  the 
corresponding  in-phase  part  is  of  such  magnitude  that  Eq.(4.23)  yields  a  non-imaginary 
flutter  frequency  (£>.  The  latter  condition  is  met  when  the  rotational  axis  is  ahead  of  the 
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one  quarter  chord  line.  Hence  by  this  method  we  can  find  the  flutter  airspeed  and 
frequency. 

The  variation  of  the  real  and  imaginary  parts  of  the  aerodynamic  moment, 
based  on  Theodorsen’s  thin  airfoil  theory,  for  k  defined  as  coc/2U_,  is  shown  in  Figure  4- 
5.  For  the  definition  ok  k  according  to  Eq.(4.2),  the  imaginary  moment  becomes  zero  at  a 
reduced  frequency  of  k=0.076.  This  is  consistent  with  the  previous  discussion  which 
showed  that  torsional  flutter  occurs  only  at  low  reduced  frequencies.  For  further  details 


Figure  4-5  Variation  with  reduced  frequency  k  of  the  real  and  imaginary  parts  of  the  dimensionless 
.  aerodynamic  moment  nio  due  to  pitching  of  an  airfoil  about  its  leading  edge  in 

incompressible  flow  (for  k  defined  as  coc/2LL). 

G.  TWO  DEGREE  OF  FREEDOM  FLUTTER 
1.  Physical  Explanation 

For  the  two-degrees-of-freedom  flutter  we  will  consider  two  cases:  a)  the  torsional 
deflection  leads  the  bending  deflection  by  90  degrees  as  shown  in  Figure  4-6  whereas  b) 
the  torsional  deflection  and  bending  are  in  phase  with  each  other  as  shown  in  Figure  4-8. 
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In  both  figures  the  solid  arrows  indicate  the  airfoil  motion  whereas  the  dotted  arrows 
indicate  the  airfoil  lift. 


For  simplicity  we  will  assume  purely  steady  state  values  for  the  lift  corresponding 
to  the  instantaneous  position  of  the  airfoil. 


Figure  4-6  Two  degrees  of  freedom  airfoil  motion  -  phase  difference  90° 


From  Figure  4-7  we  can  see  that  for  the  two  first  quarters  of  the  cycle,  both  lift  L 
and  displacement  dh  are  positive  whereas  for  the  two  last  quarters  of  the  cycle,  both  lift  L 
and  displacement  dh  are  negative.  As  a  result  over  the  whole  cycle  work  is  done  on  the 
airfoil. 

When  the  torsional  deflection  is  in  phase  with  the  bending  deflection  then  from 


Figure  4-8  Two  degrees  of  freedom  airfoil  motion  -  phase  difference  0° 


Figure  4-9,  we  can  see  that  for  the  first  quarter  of  the  cycle,  the  lift  L  and  the 
displacement  dh  are  positive  so  positive  work  is  done  on  the  airfoil.  For  the  second 
quarter  of  the  cycle,  the  lift  L  is  positive  whereas  the  displacement  dh  is  negative  so 
negative  work  is  done  on  the  airfoil.  Similarly  positive  and  negative  work  are  done  on  the 
airfoil  during  the  third  and  fourth  quarter  of  the  cycle.  As  a  result  zero  work  is  done  on 
the  airfoil. 

2.  Equations  of  motion 

The  equations  of  motion  for  the  two-degree-of-freedom  system  with  no 
mechanical  damping  shown,  for  example,  by  Fung  [Ref.  1]  pages  210-211,  can  be  written 
as: 

mh  +  Saa + hKh  =  -L  (4.25) 

and  Sah  +  Iad(  +  ocKa  =Ma  (4.26) 

where  Sa  is  the  static  moment  about  the  elastic  axis. 

7ais  the  airfoil  moment  of  inertia 

L  is  the  airfoil  aerodynamic  lift 
and  M a  is  the  aerodynamic  moment  about  the  elastic  axis 

a.  Time  domain  analysis 

For  the  incompressible  flow  case  the  equations  of  motion  (4.25)  and  (4.26) 
can  be  non-dimensionalized  using  the  same  non-dimensional  coefficients  that  were  used 
for  the  one  degree  of  freedom  case  Eq.(4.15)-(4.17). 

Then  the  system  of  equations  (4.25)  and  (4.26)  becomes: 
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and 


nth*  +  Saa ”  +  mk\h  =  -2c l  Ik 
Sah"  +  ia(Xr  +  iakl  =  2cm  / K 


where  the  primes  refer  to  differentiation  with  respect  to  nondimensional  time  x. 
.  In  a  matrix  form  the  above  system  of  equations  becomes: 


with 


[MfX}  +[kfx}={F} 


(4.27) 


For  the  compressible  flow  case  (N-S  Code)  the  equations  of  motion  (4.25) 
and  (4.26)  can  be  non-dimensionalized  using  the  free-stream  speed  of  sound,  ,  instead 
of  the  free-stream  velocity  f/„  such  that  the  dimensionless  time  can  be  defined  from 
Eq.(4.20). 

Then  the  system  of  equations  have  the  same  form  as  Eq.(4.27)  with  the 
difference  that  the  vector  {p}  now  is  defined  as 

Eq.(4.34)  can  be  rewritten  as: 

{x  }  =  [M  ]-'  y-  [M  ]-■  \k  Kx  }  (4.28) 
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Defining  {X  }  =  {f },  we  can  solve  Eq.(4.28)  as  a  system  of  two  first-order 
differential  equations 

{k}  =  [M  ]■'  {f]-  [M  ]•'  MX  } 

As  we  mentioned  for  the  one  degree  of  freedom  case,  the  integration  of  this  system  of 
equations  is  done  with  the  fourth-order  Runge-Kutta  scheme  in  the  panel  code,  and  with 
a  first-order  Euler  scheme  in  the  Euler/N-S  codes. 

b.  Frequency  domain  analysis 

Assuming  constant  amplitude  time  harmonic  oscillation  including  plunge 

and  torsion 

h  =  h0eieo' 
a  =  a  0e imt 

the  equations  of  motion  become: 

eim[-(D2Mh0-co2a0Sa+h0Kh]=L  (4.29) 

eim\ra2lJ-aQ-co2h0Sa+a0Ka\=Ma  (4.30) 

The  aerodynamic  lift  L  and  moment  M0  are  assumed  to  be  linear  functions  of  h0  and  Oo  so 
Eqs.  (4.30)  and  (4.31)  constitute  a  system  of  homogeneous  equations.  Therefore  the 
determinant  of  these  equations  must  be  set  to  zero  in  order  to  obtain  a  solution.  As 
discussed  in  [Ref.  1  and  6]  the  solution  of  the  flutter  determinant  normally  requires  an 
iterative  procedure.  For  further  details  refer  to  these  text  books. 
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H.  NS  -  UPOT  REDUCED  FREQUENCY  COMPATIBILITY 


As  mentioned  before,  the  reduced  frequencies  calculated  with  the  Euler  and  N.S 
codes  need  to  be  redefined  in  order  to  be  compatible  with  the  conventional  definition  of 
k.  This  is  done  as  follows: 

The  non-dimensional  time  used  in  UPOT  is  defined  as 

tUm 

7 upoT  ~  (4.31) 

c 

whereas  in  the  N.S  code  it  is  defined  as 


(4.32) 


c 

In  order  to  have  comparable  results  from  the  compressible  and  incompressible  flow 
solutions  we  have  to  adjust  the  reduced  frequency  values  of  the  NS  code  to  those  of  the 
UPOT  code.  From  the  definition  of  the  reduced  frequency: 


(4,33) 


we  get 


k  — 


2  it  c 

Tu~ 


where  T  is  the  period  of  airfoil  oscillation. 

For  the  results  computed  with  the  UPOT  code  and  because  of  (4.1)  we  have 

2  n  c  2 n 
Tupotc  Um  Tupot 


Similarly  for  the  NS  results  from  (4.4)  and  because  of  (4.2)  we  get 


(4.34) 


(4.35) 


2n  c  _  2n  c  _2n  1 
Tnsc  Um  Tnsc  Tns  M. 


(4.36) 
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The  last  equation  shows  that  we  have  to  divide  the  reduced  frequencies  that  come  from 
the  NS  code  (2  n!  Tm )  with  the  Mach  number  ( )  in  each  case  in  order  to  find  the 
reduced  frequency  which  is  compatible  with  the  UPOT  results. 

I.  FLUTTER  PREDICTION  WITH  UPOT  CODE 

One  practical  way  to  find  the  flutter  condition  for  a  given  airfoil  with  the  UPOT 
code  is  to  see  the  resulting  diagram  of  the  moment  coefficient  with  respect  to  AOA.  The 
work  is  proportional  to  the  integrated  area  in  the  loop.  Positive  area  means  that  work  goes 
into  the  flow.  Negative  area  means  that  work  goes  into  the  airfoil.  Zero  area  means  that 
the  airfoil  flutters.  Running  the  code  for  low  values  of  reduced  frequency  it  is  seen  that 
the  aerodynamic  moment  coefficient  forms  a  clockwise  elliptic  path  as  the  AOA  changes 
during  the  airfoil  oscillation  (Figure  4-10),  so  net  positive  work  is  done  on  the  airfoil. 


Figure  4-10  Aerodynamic  moment  coefficient  graph  for  low  k  values 
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For  high  values  of  reduced  frequency  the  path  has  a  counterclockwise  direction  as 
shown  in  Figure  4-11  and  the  work  is  negative. 


Figure  4-11  Aerodynamic  moment  coefficient  graph  for  high  k  values 

At  last,  when  the  critical  value  of  reduced  frequency  is  used  the  path  becomes  a 
straight  line  as  shown  in  Figure  4-12,  so  no  work  is  done  on  the  airfoil. 


Figure  4-12  Aerodynamic  moment  coefficient  graph  for  the  critical  flutter  k  value 
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V.  RESULTS 


A.  GENERAL 

All  the  results  presented  in  this  thesis  were  obtained  with  the  time  domain  method 
described  in  Chapter  IV,  and  the  flow  field  solutions  at  each  time  step  were  computed 
based  on  the  Navier-Stokes,  Euler  or  unsteady  potential  flow  (UPOT)  solution,  described 
in  Chapter  m. 

The  reduced  frequency  values  used  in  the  Navier-Stokes  and  Euler  codes  were 
redefined  in  order  to  be  compatible  with  the  UPOT  definition.  The  grid  generation  for  the 
airfoils  and  the  input  data  for  the  Navier  -  Stokes  code  are  presented  in  Sections  B  and  C. 
The  procedure  for  extracting  the  results  is  described  briefly  in  Section  D.  Then  the 
calculation  for  the  non-viscous  (Euler)  case  is  presented  in  Section  E.  The  calculation  for 
a  wide  range  of  pivot  points  is  shown  in  Section  F.  The  effect  of  Mach  number  and  airfoil 
thickness  on  torsional  flutter  is  presented  in  Sections  G  and  H.  The  influence  of  the 
viscous  effects  and  a  comparison  between  viscous  and  non-viscous  flow  are  shown  in 
Sections  I  and  J.  In  Section  K  a  comparison  between  UPOT  and  NS  results  for  low 
subsonic  Mach  numbers  (M=0. 1-0.3)  is  presented. 

B.  GRID  GENERATION 

Airfoil  grids  for  NACA  0006,  0009,  0012  and  0015  were  computed  with 
ALGEM.  In  Figure  5-1  the  NACA  0015  grid  for  the  Euler  calculations  is  presented 
whereas  in  Figure  5-2  the  grid  details  close  to  the  airfoil  are  shown.  This  NACA  0015 
coarse  grid  size  is  201X41  points. 

For  viscous  flow  different  grids  were  used  because  of  the.  necessity  to  resolve  the 
flow  close  to  the  airfoil.  The  grid  size  is  201X61  points. 
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Figure  5-1  C-grid  for  NACA  0015  (non-viscous  flow  -  every  other  grid  line  shown) 
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Figure  5*2  NACA  0015  C-grid  details  (non*viscous  flow) 

C.  INPUT  DATA  DESCRIPTION 

The  input  data  file  includes  many  variables,  some  of  which  have  to  be  changed 
during  the  calculation  procedure.  A  typical  input  file  for  the  initial  calculation  of  the 
steady  solution  for  the  non- viscous  flow  is  presented  next.  The  use  of  the  input  file 
variables  is  also  explained  in  this  section. 


INPUT  FILE : 

#  IREAD  ITER  NPRINT  NLOAD  ODVAR 

0  10000  1  1  1.00 

#  ALPHA  OSCIL  RAMP  REDFRE  ALFAMND  ALFAMXD 

0.5  false  false  0.100  -.500  0.500 

#  PLUNGE  PLMX  PLMY  PLPHSXD  PLFREQ 

false  0.  0.10  0.  0.1 

#  MACH  RE  VISC  TURBL 

0.700  0.  false  false 

#  TIMEAC  COUR  NEWTTT 

false  30.00  1 

#  free  mass  ialpha  ka  kh  xp  xa 

false  0.0  100.0  0.200  0.0  0.0  0.0 

#  aneut  hneut  hO 

0.0  0.0  0.0 

INPUT  VARIABLES  EXPLANATION: 

IREAD  0  No  initial  solution,  free  stream  conditions  initialize  the 
flowfield  in  the  startup  steady-state  computations 
1  Initial  solution  is  read  from  a  binary  file  saved  from 

the  previous  run  (default,  at  the  end  of  each  run,  solution 
file  is  saved  as  binary) 

2  Initial  solution  is  read  as  formatted  (plot3d  form) 

-1  Initial  solution  is  read  as  binary,  unsteady  motion  starts 

-2  Initial  solution  is  read  as  formatted  (plot3d  form),  . 

unsteady  motion  starts 

ITER  #  of  timesteps 

NPRINT  Residuals  are  printed  out  at  every  nprint  timesteps 
NLOAD  Aerodynamic  loads  are  printed  out  in  ns.out  and  written  into 
lo.d  file  at  every  nload  timesteps 

ODVAR  Solution  variables,  q  array,  are  written  into  "qp.d"  file 
at  every  delta  odvar  change  in  unsteady  motions 
(degrees  in  oscillatory  motion,  amplitude  change  in  plunge) 
ALPHA  Steady  state  AOA  (do  not  set  it  to  zero,  instead  set  it  to  0.0001) 
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OSCIL 

RAMP 

REDFRE 

ALFAMND 

ALFAMXD 

PLUNGE 

PLMX 

PLMY 

PLPHSXD 

PLFREQ 

REYNOLD 

MACH 

vise 

TURBL 

TIMEACC 

COUR 


false  N/A 

true  sinusoidal  oscillations  in  pitch 
false  N/A 

true  straight  ramp  motion  in  pitch 

Reduced  frequency  of  the  unsteady  pitching  motion, 
based  on  the  half  chord,  chord  length  is  assumed  to  be  1. 

Min  AOA  of  the  pitching  motion 

Max  AOA  of  the  pitching  motion 

false  N/A 

true  plunging  motion 
Plunge  amplitude  in  x 
Plunge  amplitude  in  y 

The  phase  angle  between  x  and  y  amplitudes,  in  degrees 

The  reduced  frequency  of  the  plunging  motion, 
based  on  the  half  chord,  chord  length  is  assumed  to  be  1. 

Reynolds  number  of  the  freestream  flowfield 

Mach  number  of  the  freestream  flowfield 

false  Euler  solution 

true  Viscous  Navier-Stokes  solution 

false  Laminar  flow  is  assumed 

true  Baldwin-Lomax  turbulence  model  is  applied 

false  Variable  local  time  stepping  in  the  computational  grid,  used 
in  the  computations  of  steady  state  and/or  attached  flowfields. 
true  Constant  time  stepping  everywhere  in  the  computational  grid, 
used  in  the  computations  of  unsteady  and  separated  flowfields 

Courant  number  of  the  timestepping  (50-1500),  determines  the 
time  step  of  the  computations  based  on  the  minimum  grid  size  and 
the  freestream  conditions.  Its  value  depends  on  the  computational  grid. 
For  diverging  computations,  its  value  needs  to  be  reduced. 
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If  the  residuals  in  the  output  file  increases  in  time,  it  is  the  sign  that 
Courant  number  is  to  be  reduced. 

NEWTIT  Number  of  Newton  subiterations  in  each  timestep,  applied  in 
unsteady  flows  (2-3),  for  steady  flowfields  it  is  set  to  1. 

D.  CALCULATION  PROCEDURE 

For  every  flutter  calculation  initially  the  steady  solution  had  to  be  found.  For  the 
non-viscous  cases  the  code  was  run  for  6000  iterations  (1000  with  COUR  Number  of  1, 
1000  with  COUR  Number  of  2,  2000  with  COUR  Number  of  4  and  1000  with  COUR 
Number  of  8).  For  the  viscous  case  and  due  to  the  increased  number  of  nodes  of  the  grid 
the  calculation  needed  over  12000  iterations  in  order  to  converge  to  a  steady  solution. 

The  steady  solution  was  assumed  to  be  satisfactory  when  the  step  variation  of  Cl 
coefficient  was  found  to  be  less  than  10'5  or  when  the  five  first  digits  of  Cl  remained 
constant. 

After  the  steady  solution  was  found  the  oscillating  airfoil  calculation  started. 

From  Chapter  m  Section  G  we  saw  that  for  pivot  points  behind  the  quarter  chord  point 
(0.25c)  a  statically  stable  solution  is  always  obtained.  Therefore  the  calculations  were 
performed  for  pivot  points  ahead  of  the  quarter  chord  point.  The  pivot  points  that  were 
chosen  for  the  flutter  calculations  were:  0.15-c,  0,  -0.25-c,  -0.50  c,  -0.75-c  and  -0.85-c.  For 
e  ery  case  an  initial  value  of  the  uncoupled  reduced  natural  pitching  frequency  k«  (Eq. 
4.3)  was  assumed  (which  corresponds  in  Eq.(4.4)  to  an  airfoil  spring  constant  K«)  and  the 
time  variation  of  the  airfoil  angle  of  attack  (AO A)  was  computed.  If  the  AOA  was  found 
to  decrease,  the  assumed  value  of  kawas  adjusted  until  the  motion  started  to  diverge.  For 
every  computation  the  slope  of  the  AOA  amplitude  was  determined  and  an  interpolation 
procedure  was  implemented  in  order  to  find  the  value  of  k«  which  gives  the  critical  flutter 
condition.  After  the  critical  condition  was  found  the  reduced  frequency  of  oscillation  was 
calculated. 
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E.  NON  -  VISCOUS  (EULER)  CALCULATION 


The  whole  procedure  described  above  will  be  shown  for  the  non- viscous  flow 
over  an  oscillating  airfoil  NACA  0015  at  M=0.7.  The  elastic  axis  of  oscillation  is  the  L.E 
of  the  airfoil  (pivot  point  Xp=0).  An  initial  AOA  disturbance  of  0.5°  will  be  imposed  on 
the  airfoil  in  order  to  start  the  motion. 

The  input  file  for  the  steady  solution  is  as  follows  (4th  run  for  COUR  Number  :8): 

#  IRE  AD  ITER  NPRINT  NLOAD  ODVAR 

1  2000  1  1  1.00 

#  ALPHA  OSCIL  RAMP  REDFRE  ALFAMND  ALFAMXD 

0.5  false  false  0.100  -.500  0.500 

#  PLUNGE  PLMX  PLMY  PLPHSXD  PLFREQ 

false  0.  0.10  0.  0.1  0.0  0.0 

#  MACH  RE  VISC  TURBL 

0.700  0.  false  false 

#  TIMEAC  COUR  NEWTIT 

false  8.00  1 

#  free  mass  ialpha  ka  kh  xp  xa 

false  0.0  100.0  0.00000  0.0  0.0  0.0 

#  aneut  hneut  hO 

0.0  0.0  0.0 

The  steady  solution  is  presented  in  Figure  5-3.  It  is  seen  that  Cl  reaches  a  steady  value  of 
0.1017  after  x  =5  sec. 


The  final  result  for  the  steady  solution  is  presented  in  Table  5-1. 


X 

AOA 

Cl 

Cd 

Cm 

IBgama 

0.5 

0.1017382 

0.0040226 

-0.0257816 

Table  5-1  Steady  state  solution  for  NACA  0015  M=0.7  (Euler  case) 
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Figure  5-3  Steady  state  solution  for  NACA  0015  M=0.7  (Euler  case) 

After  the  steady  solution  is  calculated,  the  critical  value  of  the  reduced  natural 
pitching  frequency  ka  for  flutter  is  found,  by  first  assuming  an  initial  value  of  0.36,  as 
shown  in  the  input  file 


# 

IREAD 

ITER 

NPRINT 

NLOAD 

ODVAR 

-1 

10000 

1 

1 

1.00 

# 

ALPHA 

OSCIL 

RAMP 

REDFRE 

ALFAMND 

ALFAMXD 

in 

o 

false 

false 

0.100 

-.500 

0.500 

# 

PLUNGE 

PLMX 

PLMY 

PLPHSXD 

PLFREQ 

false 

0. 

o 

H 

o 

0. 

0.1 

o 

o 

# 

MACH .  ■ 

RE 

•  vise 

TURBL 

0.700 

0. 

false 

false 

# 

TIMEAC 

COUR 

NEWT  IT 

true 

15.00 

3 
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# 

free 

mass 

ialpha 

ka  kh 

xp 

xa 

true 

o 

o 

100.0 

0.360000  0.0 

o 

o 

o 

o 

# 

aneut 

hneut 

hO 

o 

o 

o 

o 

o 

o 

As  one  can  see,  seven  variables  changed  values  comparing  with  the  steady  case 
input  file:  the  IREAD  was  turned  to  -1,  the  number  of  iterations  to  10000,  the  TIMEAC  to 
true,  the  COUR  number  to  15.0 ,  the  NEWTIT  to  3  and  the  free  to  true.  Also  ialpha  (non- 
dimensional  moment  of  inertia)  was  chosen  to  be  100  and  the  ka  to  be  0.36. 

The  plot  of  the  resulting  airfoil  motion  is  presented  in  Figure  5-4  . 


From  the  plot  we  can  see  that  the  airfoil  motion  decays  with  time  so  the  assumed 
value  of  k«  was  too  high.  In  order  to  determine  the  critical  ka  value  the  maximum  AOA 
values  of  the  graph  are  input  into  a  regression  equation  to  find  the  slope  of  the  maximum 
AOA  values,  namely  -2.860E-04  (Figure  5-5): 
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Figure  5-5  Time  rate  of  change  of  AOA  amplitude  of  airfoil  motion  K„=0.36  (Euler  case) 


Repeating  this  process  for  ka  of  0.25  we  get  the  AOA  history  presented  in  Figure 


From  this  plot  we  observe  that  the  airfoil  motion  diverges.  Therefore  the  assumed 
value  of  ka  is  too  low.  Then  applying  again  a  regression  equation  for  the  maximum  AOA 
values  of  the  graph  yields  a  slope  of  4.566E-04  (  Figure  5-7): 


Figure  5-7  Time  rate  of  change  of  AOA  amplitude  of  airfoil  motion  K^O.25  (Euler  case) 

Now  we  can  interpolate  the  values  of  the  slope  found  for  the  two  assumed  values 
of  the  reduced  natural  pitching  frequency  ka.  In  order  to  have  more  accurate  results  we 
can  follow  the  same  procedure  for  another  value  of  k«  (here  assumed  0.30).  Finally  we 
get  the  results  shown  in  Table  5-2. 


Ka 

Slope  1 

0.36 

-2.860E-04 

0.30 

1.669E-04 

0.25 

4.5664E-04 

Table  5-2  Regression  equation  slope  values  for  various  ka  values 


Making  a  curvefit  in  EXCEL  for  the  above  results  (Figure  5-8)  we  can  find  that 
the  value  of  the  ka  for  which  the  slope  of  the  extreme  points  is  zero,  is  ka=0.324. 
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y  =  -0.0068X  +  0.0022 


Figure  5-8  Curvefit  of  Ka  vs  time  rate  of  change  of  AOA  amplitude  (Euler  case) 

To  verify  whether  this  ka  value  is  indeed  the  critical  reduced  natural  pitching 
frequency,  the  code  was  run  for  k«  =0.324  which  yielded  the  airfoil  motion  presented  in 
Figure  5-9: 


Figure  5-9  Unsteady  case  solution  -  Ka=0.324  Critical  flutter  case  (Euler) 
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Omitting  the  first  peaks,  which  are  caused  by  the  starting  transient,  and  applying 
again  a  regression  equation  for  the  maximum  points,  we  find  that  the  slope  is 
6.638591-10'7,  which  is  very  close  to  zero  (Figure  5-10).  Therefore  the  reduced  natural 
pitching  frequency  k«  of  0.324  is  accepted  as  the  critical  value.  If  this  value  didn’t  give 
the  flutter  condition  (slope  not  close  to  zero),  the  new  slope  value  would  be  interpolated 
again  until  finding  the  k«  for  zero  slope. 
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Figure  5-10  Time  rate  of  change  of  AOA  amplitude  of  airfoil  motion  1^=0.324  Critical  flutter 
.  (Euler  case) 

The  reduced  frequency  of  oscillation  is  found  using  (4.36) 

k  -  2 —  l 

In  order  to  find  the  period  of  oscillation  the  initial  transient  exhibiting  variable 
peaks  are  ignored,  and  the  average  period  between  the  minimum  and  maximun  peaks  is 
computed  as  shown  in  Table  5-3. 


59 


Min  peaks 

T  (NS) 

Max  peaks 

T  (NS) 

92.6997 

105.954 

119.209 

26.5092 

132.465 

26.51039 

145.721 

26.51206 

158.977 

Average  T  (NS)  = 

26.511 

Table  5-3  Average  period  T  calculation  for  NACA  0015  flutter  -  Ka=0.324 


The  average  period  of  oscillation  (non-dimensionalized)  is  26.51 1.  Then  applying 
Eq.(4.36)  for  Mach  number  of  0.7  we  find  that  the  reduced  frequency  of  oscillation  is 
0.339. 


The  pressure  distribution  around  the  airfoil  for  a  whole  cycle  of  oscillation  at  30° 
increments  is  shown  in  Figure  5-11. 
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AOA=150‘ 


AOA=180° 

Figure  5-11  Pressure  Contours  around  NACA  0015  for  a  whole  cycle  of  oscillation  M=0.7  Xp=0.0 
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AOA=330‘ 


AOA=360' 


Fig.5-11  Pressure  Contours  around  NACA  0015  for  a  whole  cycle  of  oscillation  M=0.7  Xp=0.0  (cont) 
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F.  EFFECT  OF  PIVOT  POINT  ON  TORSIONAL  FLUTTER 


The  previously  described  procedure  was  applied  to  obtain  the  effect  of  pivot 
location,  Mach  number  and  airfoil  thickness  on  single  degree  of  freedom  torsional  flutter. 
Figure  5-12  shows  the  variation  of  the  flutter  reduced  frequency  with  pivot  point  for  the 
NACA  0015  airfoil  at  M=0.7. 


Figure  5-12  Effect  of  Pivot  Point  on  reduced  frequency  for  NACA  0015  and  M=0.7  (Euler  case) 


The  area  below  the  curve  is  an  unstable  region  for  the  airfoil;  that  is  if  the 
frequency  of  oscillation  is  lower  than  the  value  corresponding  to  the  curve  for  a  certain 
pivot  point,  then  flutter  is  triggered  by  an  initial  disturbance  of  0.5  degrees.  In  contrast 
the  area  above  the  curve  is  the  stable  region  and  every  oscillation  will  die  out. 

It  is  readily  seen  that  the  most  flutter  prone  range  of  pivot  points  is  between  0  and 
-0.50  c.  The  pivot  point  most  likely  to  induce  flutter  is  Xp=-0.25-c,  i.e  when  the  pivot 
point  of  oscillation  is  one  quarter  chord  ahead  of  the  L.E. 

The  above  information  can  be  shown  in  terms  of  the  reduced  natural  pitching 
frequency  k«  required  to  avoid  flutter.  It  should  be  noted  that  for  elastic  axis  locations  of 
-0.85<Xp<0.15  the  used  value  of  moment  of  inertia  iawas  100.  For  Xp=0.15  and  -0.85 


63 


the  solutions  were  always  decaying  with  this  value  of  i«  so  an  increased  value  of  300  was 
used.The  resulting  graph  is  shown  in  Figure  5-13: 
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Figure  5-13  Effect  of  Pivot  Point  on  critical  kot  for  NACA  0015,  M=0.7  and  ia=100  (Euler  case) 


Again,  the  area  below  the  curve  is  an  unstable  region  for  the  airfoil;  if  the  airfoil 
spring  constant  K«  is  lower  than  the  value  that  corresponds  to  the  curve  for  a  certain  pivot 
point  (low  torsional  stiffness),  then  flutter  is  induced  by  an  initial  disturbance  of  0.5  rees. 

G.  EFFECT  OF  MACH  NUMBER  ON  TORSIONAL  FLUTTER 

Figure  5-14  expands  the  information  presented  in  Figure  5-12  by  displaying  the 
effect  of  Mach  number  on  torsional  flutter. 

An  increase  in  Mach  number  significantly  increases  the  region  of  instability  for 
the  whole  range  of  pivot  points. 

Another  important  result  that  comes  from  this  graph  is  that  the  values  of  reduced 
frequency  tend  to  become  constant  when  reducing  the  Mach  number  to  very  low  subsonic 
speeds  of  M=0. 1-0.2. 
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NACA  0015 


Xp 

Figure  5-14  Effect  of  Mach  Number  and  Pivot  Point  for  NACA  0015  (Euler  case) 


This  information  is  shown  in  bar  chart  form  in  Figure  5-15: 


Figure  5-15  Reduced  frequency  variation  with  Mach  number  for  Pivot  Points  -0.85<Xp<0.15 
(NACA  0015) 
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The  percentage  increase  of  the  reduced  flutter  frequency  with  Mach  number  and 
pivot  points  is  tabulated  in  Table  5-4. 


M  \  Xp 

0.15 

0 

-0.25 

-0.5 

-0.75 

-0.85  | 

0.7 

12.4% 

34.6% 

59.8% 

68.6% 

45.8% 

0.5 

11.4% 

13.8% 

23.4% 

28.9% 

28.3% 

0.3 

7.4% 

1.9% 

4.9% 

6.9% 

13.3% 

0.2 

2.4% 

-1.1% 

-2.1% 

-2.0% 

-0.1% 

4.0% 

0.1 

0 

0 

0 

0 

0 

0 

TaWe  5-4  Percentage  Increase  of  k  with  Mach  number  and  Pivot  Points 


Another  way  to  display  the  effect  of  pivot  location  and  Mach  number  is  to  plot  the 
reduced  natural  torsional  frequency  ka  needed  to  suppress  flutter.  Again  it  is  seen  from 
Figure  5-16  that  an  increase  in  Mach  number  is  destabilizing. 


Xp 


Figure  5-16  Effect  of  Mach  number  and  pitch  axis  location  on  torsional  flutter  of  NACA  0015 
airfoil,  i„=100 
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Another  useful  way  of  representing  the  results  is  to  show  the  variation  of  the 
reduced  flutter  velocity  with  Mach  number  for  every  pivot  point.  As  discussed  in  Chapter 

HI  the  reduced  flutter  velocity  is  defined  as  V*.  =  =  —  and  can  be  used  in  order  to 

find  the  speed  U  for  which  flutter  occurs  for  given  values  of  K«.  Figure  5-17  shows  the 
decrease  of  the  flutter  speed  as  the  Mach  number  is  increased.  The  results  shown  are  only 
for  pivot  points  -0.85<Xp<0  for  which  ia=100. 


Xp 


Figure  5-17  Variation  of  flutter  speed  with  Mach  Number  for  NACA  0015,  io=100 

H.  EFFECT  OF  AIRFOIL  THICKNESS 

The  same  procedure  was  applied  to  study  the  effect  of  airfoil  thickness. 
Calculations  were  done  for  the  airfoils:  NACA  0006,  NACA  0009,  NACA  0012  and 
NACA  0015  for  a  Mach  number  of  0.7.  The  variation  of  critical  reduced  frequency  with 
the  airfoil  thickness  is  presented  in  Figure  5-18  and  Figure  5-19  showing  that  over  the 
whole  pivot  point  range  any  increase  in  airfoil  thickness  has  a  destabilizing  effect. 
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Figure  5-18  Effect  of  airfoil  thickness  on  critical  reduced  frequency  at  M=0.7 


Figure  5-19  Effect  of  airfoil  thickness  at  M=0.7 
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This  destabilizing  influence  of  airfoil  thickness  is  quantified  in  Table  5-5  where 
the  percentage  increase  of  the  reduced  flutter  frequency  k  is  tabulated. 


Airfoil 

thickness 

0.15 

0 

-0.25 

-6.5 

-0.75 

-0.85 

6 

0 

0 

0 

0 

0 

0 

9 

5.6% 

6.8% 

6.3% 

mm 

3.2% 

12 

11.2% 

37.8% 

49.8% 

18.8% 

5.2% 

5.4% 

15 

28.1% 

61.8% 

70.4% 

58.6% 

14.3% 

10.3% 

Table  5-5  Percentage  increase  of  k  with  airfoil  thickness  and  pivot  point  location 


Similarly,  Figure  5-20  displays  the  reduced  natural  pitching  frequency  ka  needed 
to  suppress  flutter.  Again  it  is  seen  that  the  NACA  0015  airfoil  requires  a  much  higher 
torsional  frequency  than  the  NACA  0006  airfoil.  Again  the  results  shown  are  for  elastic 


Xp 


Figure  5-20  Effect  of  airfoil  thickness  on  torsional  flutter  atM=0.7,  i„=100 
axis  locations  of  -0.85<Xp<0.15  for  which  the  used  value  of  moment  of  inertia  was  100. 
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Finally  the  destabilizing  effect  is  also  shown  in  Figure  5-21.  Increasing  the  airfoil 
thickness  results  in  a  significant  decrease  in  the  flutter  speed. 


Xp 


Figure  5-21  Variation  of  flutter  speed  with  airfoil  thickness  at  M=0.7,  io=100 


I.  VISCOUS  (NA  VIER-STOKES)  CASE  CALCULA TION 


The  same  procedure  previously  described  for  Euler  case  will  be  shown  for  the 
viscous  flow  over  the  same  oscillating  airfoil  (NACA  0015),  the  same  Mach  number 
(M=0.7)  and  the  same  elastic  axis  of  oscillation  (L.E  of  the  airfoil  -  pivot  point  Xp=0). 
The  input  file  for  the  steady  solution  is  as  follows: 


# 

IREAD 

ITER 

NPRINT 

NLOAD 

ODVAR 

0 

20000 

1 

1 

O 

O 

# 

.  ALPHA 

OSCIL 

RAMP 

REDFRE 

ALFAMND 

ALFAMXD 

0.5 

false 

false 

0.100 

-.500' 

0.500 

# 

PLUNGE 

PLMX 

PLMY 

PLPHSXD 

PLFREQ 

false 

0. 

0.10 

0. 

0.1 

o 

o 

# 

MACH 

RE 

vise 

TURBL 

0.700 

1 . 0e6 

true 

true 
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# 

TIME AC 

COUR 

NEWTIT 

true 

1000.0 

1 

# 

free 

mass 

ialpha 

ka 

kh 

xp 

xa 

false 

0.0 

100.0 

0.00000 

0.0 

o 

o 

0.0 

#  ■ 

aneut 

hneut 

hO 

0.0 

0.0 

0.0 

The  steady  solution  is  presented  in  Figure  5-22.  We  can  see  that  Cl  reaches  a 
steady  value  of  0.08977  after  t  =40  sec. 
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Figure  5-22  Steady  state  soln  for  NACA  0015  M=0.7  (N-S  case) 

The  final  result  for  the  steady  solution  is  shown  at 
Table  5-6: 


T 

Aoa 

Cl 

Cd 

Cm 

0.5 

0.089777 

0.006215 

0.000665 

Table  5-6  Steady  state  solution  for  NACA  0015  M=0.7  (N-S  case) 
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After  the  steady  solution  is  obtained,  the  flutter  condition  can  be  found  using  the 
same  procedure  described  for  the  Euler  calculations.  We  start  with  a  pivot  point  of 
Xp=0.0  and  an  assumed  value  of  0.36  for  k«.  The  following  is  the  input  file  for  the  N.S. 
code: 


# 

IREAD 

ITER 

NPRINT 

NLOAD 

ODVAR 

-1 

30000 

1 

1 

1.00 

# 

ALPHA 

OSCIL 

RAMP 

REDFRE 

ALFAMND 

ALFAMXD 

0.5 

false 

false 

0.100 

-.500 

0.500 

# 

PLUNGE 

PLMX 

PLMY 

PLPHSXD 

PLFREQ 

false 

0. 

0.10 

0. 

0.1 

0.0  0.0 

# 

MACH 

RE 

vise 

TURBL 

0.700 

1 . 0e6 

true 

true 

#  . 

TIMEAC 

COUR 

NEWTIT 

true 

700.00 

3 

# 

free 

mass 

ialpha 

ka 

kh 

xp  xa 

true 

0.0 

100.0 

0.340000 

0.0 

0.0  0.0 

# 

aneut 

hneut 

hO 

0.0 

0.0 

0.0 

The  plot  of  the  resulting  airfoil  motion  is  presented  in  Figure  5-23. 


Figure  5-23  Unsteady  case  solution  -  Ka=0.34  (N-S  case) 
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Since  the  airfoil  motion  is  seen  to  decay  the  assumed  value  of  ka  is  too  high.  After 
finding  the  maximum  points  of  the  graph  and  applying  a  regression  equation  we  find  that 
the  slope  of  the  maximum  points  is  -3.188E-04  (Figure  5-24). 


Figure  5-24  Time  rate  of  change  of  AOA  amplitude  for  Ka=0.34  (N-S  case) 
Similarly  for  ko=0.20  the  diverging  motion  presented  in  Figure  5-25  is  obtained. 


Then  applying  a  regression  equation  for  the  maximum  points  of  the  graph,  the 
slope  is  found  to  be  4.566E-04  (Figure  5-26): 


y  =  4.383E-04X  +  5.744E-01 


Figure  5-26  Time  rate  r  hange  of  AOA  amplitude  for  Ka=0.20  (N-S  case) 


The  results  for  the  two  trials  are  shown  in  Table  5-7  which  make  it  possible  to  interpolate 
the  value  of  k«  which  produces  a  constant  amplitude  of  oscillation.  As  shown  in  Figure  5- 
27  this  value  is  k«=0.252. 


Ka 

Slope 

0.34 

-2.860E-04 

0.20 

4.383E-04 

Table  5-7  Regression  equation  slope  values  for  various  ka  values 
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y  =  -5.174236E-03X  +  1 .473201 E-03 
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Figure  5-27  Curvefit  of  Ka  vs  time  rate  of  change  of  AOA  amplitude  (N-S  case) 


Using  this  ka  value  indeed  produces  the  oscillation  shown  in  Figure  5-28  and  the 
slope  -6.783E- 10'5  shown  in  Figure  5-29. 
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Figure  5-28  Unsteady  case  solution  -  Ka=0.252  Critical  flutter  case  (N-S) 
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y  =  -6.783E-05X  +  5.466E-01 


time  (NS) 


Figure  5-29  Time  rate  of  change  of  AOA  amplitude  for  Ka=0.324  Critical  flutter  (N-S  case) 

From  the  airfoil  motion  we  find  that  the  average  period  of  oscillation  is  33.215 
(non-dimensional  time)  and  reduced  frequency  of  oscillation  0.270. 


J.  EULER-NS  RESULTS  COMPARISON 

Due  to  time  constraints  the  N-S  study  was  performed  only  for  NACA  0015  airfoil 
and  M=0.7.  The  critical  reduced  frequency  values  for  pivot  points  between 
-0.85<Xp<0. 15  are  presented  in  Figure  5-30  and  compared  with  the  Euler  predictions. 

It  is  seen  that  the  values  of  critical  reduced  frequencies  found  with  the  Euler 
calculations  are  up  to  23%  higher  than  the  N.S  predictions.  As  expected,  viscous  flow 
effects  reduce  the  possibility  of  flutter.  From  Eq.(4.5)  it  is  seen  that  the  same  airfoil  at  the 
same  Mach  number  with  60%  less  torsional  rigidity  will  have  the  same  flutter  stability. 
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Figure  5-30  Euler  -  NS  reduced  frequency  results  for  NACA  0015  and  M=0.7 


Table  5-8  shows  the  numerical  differences  between  the  above  calculations  for 
Euler  and  Navier-Stokes  cases. 


Xp 

0.15 

0.00 

-0.25 

-0.50 

-0.75 

-0.85 

k(Euler) 

0.223 

0.339 

0.380 

0.334 

0.233 

0.221 

k  (N-S) 

0.191 

0.270 

0.310 

0.260 

0.217 

0:205 

Difference 

14.3% 

20.3% 

18.4% 

22.1% 

6.8% 

7.2% 

Table  5-8  Numerical  differences  of  reduced  frequency  between  Euler  and  N-S  calculations  for 
NACA  0015  and  M=0.7 


The  reduced  natural  pitching  frequency  ka  values  computed  with  the  N-S  and 
Euler  codes  are  presented  in  Figure  5-31  indicating  that  the  minimum  torsional  stiffness 
necessary  to  prevent  flutter  is  smaller  using  viscous  flow  calculations. 
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In  order  to  explain  this  finding  it  is  instructive  to  examine  the  pressure  distri¬ 
butions  and  Mach  number  contour  plots  for  both  cases  in  steady  AOA  of  0.5°  (Figure  5- 
32  to  Figure  5-39).  The  main  difference  of  the  graphs  has  to  do  with  the  strength  of  the 
shock  above  the  airfoil.  For  the  Euler  case  it  is  seen  that  a  relatively  strong  shock  is 
formed  over  the  airfoil.  The  graphs  that  come  from  the  Navier-Stokes  calculations  show 
that  the  shock  over  the  airfoil  is  very  weak.  This  is  due  to  the  boundary  layer  which 
smooths  the  shock. 


Xp 


Figure  5-31  Euler  -  NS  reduced  natural  pitching  frequency  results  for  NACA  0015  and  M=0.7 
i„=100 
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Figure  5-32  Pressure  distribution  contours  for  NACA  0015  M=0.7  steady  AOA  (Euler) 


Figure  5-33  Pressure  distribution  contours  for  NACA  0015  M=0.7  steady  AOA  (Euler-detail) 
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Figure  5-34  Mach  contours  for  NACA  0015  M=0.7  steady  AOA  (Euler) 


Figure  5-35  Mach  contours  for  NACA  0015  M=0.7  steady  AOA  (Euler-detail) 


Figure  5-37  Pressure  distribution  contours  for  NACA  0015  M=0.7  steady  AOA  (N-S-detail) 
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Figure  5-38  Mach  contours  for  NACA  0015  M=0.7  steady  AOA  (N-S) 


82 


K.  PANEL  -  EULER  CODE  COMPARISON 


Another  interesting  area  of  study  is  the  comparison  of  results  from  the  Euler  code 
for  non-viscous  low  subsonic  flow  with  those  from  the  UPOT  panel  code  which 
calculates  inviscid  incompressible  flow  over  the  airfoil.  The  calculations  were  again 
made  for  flow  the  NACA  0015.  The  results  for  M=0.3, 0.2, 0.1  and  incompressible  flow 
are  presented  in  Figure  5-40. 


NACA  0015 


Figure  5-40  Low  subsonic  Euler  code  results  comparison  with  UPOT  results. 

It  is  seen  that  the  trends  predicted  by  both  codes  are  in  good  agreement,  but  there 
is  a  significant  quantitative  difference  between  the  Euler  prediction  at  M=0.1  and  the 
incompressible  panel  code  prediction.  This  difference  could  be  reduced  substantially  if 
the  panel  code  was  run  with  a  very  small  time  step,  as  shown  in  Figure  5-40  and  Figure 
5-41. 
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%  difference 


Figure  5-41  %  differences  between  Euler  and  UPOT  results  (w.r.t  M=0.1  results) 


8. 


VI.  SUMMARY 


The  time  domain  flutter  analysis  of  NACA  0006,  NACA  0009,  NACA  0012  and 
NACA  0015  airfoils  presented  in  this  thesis  leads  to  the  following  major  conclusions: 

A.  Linearized  incompressible  and  linearized  subsonic  compressible  flow 
theory  yields  unconservative  estimates  of  single-degree-of-freedom  torsional  airfoil 
flutter. 

B.  Torsional  flutter  region  increases  with  increasing  airfoil  thickness  and 
Mach  number. 

C.  Viscous  flow  effects  have  a  stabilizing  influence. 

D.  Pitch  axis  locations  upstream  of  the  quarter  chord  point  may  induce  flutter, 
especially  in  the  range  between  the  leading  edge  and  a  point  a  half-chord  upstream  of  the 
leading  edge. 

E.  Flutter  boundaries  computed  with  the  incompressible  panel  code  UPOT 
require  very  small  time  steps;  however,  even  for  the  smallest  time  step  used,  the  UPOT 
computed  flutter  boundaries  differed  from  Euler  computed  boundaries  by  a  significant 
amount. 

F.  Euler  computations  converged  relatively  quickly  for  high  subsonic  Mach 
numbers,  but  required  increasingly  large  times  for  completion  as  the  Mach  number  was 
•decreased  toward  0.1.  Every  Euler  computation  could  be  completed  in  about  12  to  24 
hours. 

G.  Navier-Stokes  computation  times  typically  were  4  to  6  times  longer  on  the 
Department’s  SGI  work  stations. 
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VII.  RECOMMENDATIONS 


It  will  be  interesting  to  study  the  effect  of  the  following  parameters: 

A.  Reynolds  Number 

Extend  the  calculation  to  Reynolds  numbers  other  than  1-106 

B.  Mach  Number 

Extend  the  calculations  to  cover  the  transonic  Mach  number  range  and 
investigate  the  effect  of  shock  motion  and  shock  boundary  layer  interaction  on 
flutter 

C.  Moment  of  Inertia 

Investigate  the  effect  of  moment  of  inertia  variation  on  flutter 

D.  Two-Degree-of-Freedom  Flutter 

Extend  the  single-degree-of-freedom  analysis  to  the  bending-torsion 
flutter  problem  and  investigate  the  effect  of  airfoil  thickness  on  this  type  of  flutter. 
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